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TECHNICAL MEMORANDUM 


Mars Global Reference Atmospheric Model (Mars-GRAM Version 3.34): Programmer' 

Guide 


1. INTRODUCTION 


1.1 Background 

The Mars Global Reference Atmospheric Model (Mars-GRAM) was developed 
(Johnson et al., 1989; Justus 1990, 1991; Appendix B, version 2.21) as an engineering- 
oriented, empirical model of the Mars atmosphere. The model is based on surface and 
atmospheric temperature data observed during the Mariner and Viking (orbiter and lander) 
missions and on surface pressure data observed by Viking landers. At higher altitudes 
(above about 120 km), Mars-GRAM is based on the Stewart (1987) thermospheric model 
given in Pitts et al., (1990?). The model provides both mean and mountain-wave 
perturbed atmospheric density for any location (height, latitude, longitude) and time 
(seasonal, diurnal). Other atmospheric variables include atmospheric temperature, 
pressure, and wind components. Dust storm effects, controlled by user-selected options, 
are provided for the atmospheric parameters. 

A second official release of Mars-GRAM (James and Justus, 1993; Appendix C, 
version 3.1) added several new capabilities, e.g., the option to simulate either local-scale 
or global-scale dust storms and density perturbations from tidal waves by the Zurek wave 
model given in Pitts et al. (1990?). 

Table 1-1 gives a brief history of Mars-GRAM program development through the 
current version — 3.34. Newest features include (1) a limitation based on atmospheric 
stability considerations for magnitude of mountain-wave density perturbations, (2) 
comparisons of density, temperature, and pressure with the COSPAR reference 
atmosphere (Pitts et al., 1990?), and (3) a new method for estimating the diurnal range of 
surface temperature based on diurnal variability of surface-absorbed solar energy. 
Technical descriptions of these new features are discussed in Section 3.0. 



Table 1-1. History of Mars-GRAM Program Versions (Page 1 of 2) 


Version 

1.00 

2.00 

2.10 

2.11 

2.20 

2.21 

2.22 

3.0 

3.1 


3.1 

3.2 


Date 


Comments 


5/20/88 

7/1/89 

9/2/89 

10/2/89 

10/7/89 

10/8/89 

11/16/89 

10/14/91 


12/17/92 


3/14/94 

11/28/94 


Preliminary version with earlier Stewart thermosphere 
model and no realistic latitude- longitude variation. 
Documented in ED44-5-20-88 preliminary report. 

Version documented in July, 1989 technical report. Has 
newer Stewart thermosphere model and realistic latitude- 
longitude variability . 

Adds version numbers to main and all subroutines. 

Corrects formats 790 and 795 in DATASTEP subroutine. 

Corrects "ATIO" to -RATIO" and "FH" to " PFH" in THERMOS 
subroutine. Corrects MARSGRAM, ATMOS 2 , DATASTEP, 

PRESSURE, P SURF ACE, STEWART 2 , STRATOS, TEMPS, and 
TSURFACE to have lines <= 72 characters in length. 

Corrects illegal log call in THERMOS by adding ES 
factors to ZF in ATMOS 2 . Adds EScalc subroutine in 
STEWART 2 . Changes name of terrain height file to 
HEIGHTS . DAT 

Removes character data from COMMON DATA and puts it in 
new COMMON FILENAME. 

Adds REAL J2 to RELLIPS, and REAL nmals to PSURFACE. 
Changes Julian date by -0.5 to be consistent with 
astronomical convention of day starting at Greenwich 
noon. ORBIT adds back 0.5 to Julian day for consistency 
with derivation of coefficients. 

Adds option for local-scale dust storm, and Zurek wave 
perturbation model. Allows heights to go "below" local 
terrain height and return "realistic" pressure, density 
and temperature, not the surface values. Both 
Interactive and Batch versions available. Batch version 
uses NAMELIST input, and is completely modular, so main 
driver program can easily be replaced by any calling 
program, such as a trajectory simulation program. 

Change comments and code for DENSRP output to 
perturbations in %; change DENSLO output file to DENSRM, 
containing random perturbation magnitude in %; change 
DENSHI output file to DENSWA, containing wave 
perturbation in %; add DENSWA to the OUTPUT file. 

Modify DENSHI and DENSLO to include the wave 
perturbation amplitude. Delete several unused variables 
from declaration statements in MAIN routine and 
subroutines SETUP, ATMOS 2 , PSURFACE, STEWART 2 , DZDUST, 
THERMOS , and STRATOS . 

Transferred version 3.1 to MSFC Unix environment and 
tested . 

Corrected DENSHI, DENSLO initial value problem. Added 
parameter value output to DATASTEP. Modified DZDUST to 
use same dust storm start time and intensity as lower 
altitude. Changed output file units to iuO for messages 
(normally screen) and iup for LIST output (iup = 0 also 
suppresses LIST and other output in batch version) . 
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Table 1-1. History of Mars-GRAM Program Versions (Page 2 of 2) 

Version Date Comments 


Corrected error in calculation of time-step 
correlations . Added line number codes in columns 73-80. 
Sorted MARSSUBS into alphabetical order by line number 
codes. Renamed Commons RAND to RANDCOM and DATA to 
DATACOM. Modified E- format in LIST output to have a 
leading digit. Made code consistent for indentation of 
If .. Then .. Else segments and Do.. Continue loops. Moved 
NAMELIST read into SETUP for batch version and removed 
commons from the batch main {to simplify use of SETUP 
and DATASTEP as subroutines in other driver programs) . 
Deleted additional unused variables in ATMOS2 and ORBIT. 
Defined pil80 in DATASTEP and corrected NVAR to NVARX. 
Corrected slight inconsistency in pressure interpolation 
in subroutine STRATOS . Added time (rel. to initial 
time) in sols to standard LIST file. Added time and 
solar longitude to OUTPUT list. 

3.3 2/7/95 Simplified vertical interpolation method in subroutine 

STRATOS. Added DENSLO and DENSHI output files back ("1 
sigma" density envelopes) . Corrected erroneous comments 
associated with file open statements in batch version 
code. Added maxfiles option to suppress output of TMAX, 
TMIN and TAVG files for systems that cannot have more 
than 16 files open at one time (i.e. set maxfiles = 16 
in Block Data routine) . 

3.31 3/28/95 Added check for density perturbations not to exceed 

value set by wave instability (i.e. maximum perturbation 
magnitude consistent with d-theta/dz > 0) . Ensured wave 
amplitudes are zero if wave perturbation model is not 
selected. Added output of F10.7 at 1AU and Mars orbit. 

3.32 4/11/95 Added COSPAR NH mean atmosphere and output of density 

deviations from COSPAR values (also new option logscale 
=2). Suppress output of NSWIND and EWWIND if maxfiles 
= 16. 

3.33 8/9/95 Takes polar cap radius from Alb function and makes it a 

separate function, polecap. Uses new regressions for 
Tavg versus Absorb (the daily average surface-absorbed 
solar flux) . Uses new regression for diurnal 
temperature amplitude (Tamp) versus Qa [the diurnal 
range (midnight-noon difference) in surface-absorbed 
solar flux] . These changes in temperature regressions 
are all in the Tsurface subroutine. 

3.34 11/1/95 Added comments for JPL Programmer’s Guide. Corrected 

logic error in density to OUTPUT file (added format 
810) . Corrected problem with crlat and polar in 
Tsurface. Corrected semi-diurnal term in dusty-case wave 
model . Corrected term in maximum mountain-wave 
perturbation model. 
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2. OVERVIEW OF THE MARS-GRAM PROGRAM 


The Mars-GRAM program consists of four FORTRAN source code files: 

(1) The Interactive form main program is MARSGRAM.FOR (or marsgram.f on UNIX 
machines) 

(2) The Batch (subroutine) form is MARSGRMB.FOR (or marsgrmb.f) 

(3) The Setup subroutine, used only in the batch form, is SETUP.FOR (or setup.f) 

(4) All other subroutines are MARSSUBS.FOR (or marssubs.f). 

Throughout this report, names of programs, subroutines, and variables appear either as all- 
capitals or as upper/lower case names. The FORTRAN compiler option used is assumed 
to be case-insensitive. 

Since version 3.2, all program code now has line numbers in columns 73-80. Line 
numbers consist of a four-character code (MARS for the interactive main, MARB for the 
batch main, SETU for the setup subroutine, etc.) and a three-digit number. For code 
inserted since version 3.2, an additional single letter appears in column 80. 

Table 2-1 gives a simplified "map" of the Mars-GRAM programs by indicating 
names of subroutines and functions, 4-character code used on source file line numbers, 
and calling subroutine(s). Appendix A gives complete descriptions of the interactive and 
batch main programs, setup subroutine, and other subroutines. Except for the setup 
subroutine, all functions and subroutines are listed in Appendix A in alphabetic order by 4- 
character code. 

Many program variables are passed between routines via "common" blocks. Table 
2-2 lists the seven common blocks and subroutines that use them. Descriptions in 
Appendix A include definitions of all input, output and local variables, except those pass 
through these common blocks. Table 2-3 lists variables in the common blocks, including 
common block name, variable name, variable type, and description of the variable. 

The interactive form of the program prompts the user for all input (with exception 
of data files HEIGHTS.DAT and COSPAR.DAT). The batch form reads and loads all 
data by the setup subroutine. The combination of setup routine and subroutines in 
marssubs.f are particularly suitable for adaptation as subroutines in user programs (e.g., 
orbital propagator and trajectory dynamics programs). See discussion of the batch form 
(MARB) in Appendix A and in Section 4. 



Both the interactive and batch forms are essentially drivers for the DATASTEP 
(and in turn the ATMOS2) subroutine (details in Appendix A). A brief outline of the 
process steps in Mars-GRAM subroutines follows: 

DATASTEP 

1 . Evaluate terrain height via TERRAIN subroutine 

2. Obtain orbital and solar positions via ORBIT subroutine 

3. Compute dust storm factors via DUSTFACT subroutine 

4. Obtain pressure gradients for winds: call ATMOS2 at the current location ± 2.5° 
latitude and longitude 

ATMOS2 

5. Calculate surface temperature via TSURFACE subroutine 

6. Compute lapse rates of atmospheric temperature via GAMMA subroutine 

7. Evaluate temperatures at "significant levels" (0, 5, 15, 30, 50 and 75 km) via TEMPS 
subroutine 

8. Obtain surface pressure via PSURFACE routine 

9. Calculate atmospheric pressure at significant levels via PRESSURE routine 
i.e., if desired height is: 

10. Below 75 km, interpolate between significant levels 

1 1. Above the base of the thermosphere, use STEWART2 

12. Between 75 km and the base of the thermosphere, interpolate via STRATOS 
subroutine 

DATASTEP 

13. Use pressure gradients, density, and Coriolis parameter to obtain wind components 
by the areostrophic (thermal wind) relations 

14. Apply viscous damping correction to areostrophic wind components 
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15. Compute random (mountain wave) perturbations in density (maintain proper 
correlation with previous perturbation value) 

16. Verify that random perturbation does not exceed stability limit 

17. Compute Zurek (tidal) wave model density perturbation 

18. Add perturbations to mean value of selected perturbation model. 
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Table 2-1. Map of Mars-GRAM Programs and Names of Subroutines, Functions, 4 

character Code, and Calling Subroutines 


Subroutine Name 

Code 

Called By 

ATMOS 2 

ATM 2 

Datastep 

Block Data 

BLKD 

N/A 

cospar 

COSP 

Datastep 

Dust fact 

DSTF 

Datastep 

Datastep 

DSTP 

MarsGRAM (Main) , Marsgrmb (Main) 

DZDUST 

DZDS 

ATMOS 2 , STEWART2 

EScalc 

ESCL 

ATMOS 2 , ST EWART 2 

gamma 

GAMA 

ATMOS 2 

orbit 

ORBT 

MarsGRAM (Main) , Datastep, SETUP 

Pressure 

PRES 

ATMOS 2 

PRSEAS 

PRSE 

ATMOS 2 , Psurface, STEWART2 

Psurf ace 

PSRF 

ATMOS 2 

RELLIPS 

RLPS 

MarsGRAM (Main) , ATMOS 2 , Datastep, 
Psurface, ST EWART 2 , THERMOS, SETUP 

SETUP 

SETU 

Marsgrmb (Main) 

Stratos 

STRA 

ATMOS 2 

STEWART2 

STW2 

ATMOS2 

THERMOS 

THRM 

ST EWART 2 

Temps 

TMPS 

ATMOS 2 

Tsurf ace 

TSRF 

ATMOS 2 

Wavepert 

WAVE 

Datastep 


Function Name 

Code 

Called By 



Alb 

ALBL 

Tsurf ace 



amp in t 

AMPN 

Wavepert 



Cp 

CPOT 

ATMOS 2 



phasint 

PHSN 

Wavepert 



polecap 

POLC 

Alb, Tsurface 



PPND 

PPND 

Datastep 



Random 

RAND 

MarsGRAM (Main) , 

Datastep , 

SETUP 

Tdiurnal 

TDIR 

Tsurface 



Terrain 

TERN 

MarsGRAM (Main) , 

Datastep, 

SETUP 
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Table 2-2. Common Blocks and Subroutines 


Common 


Used By 


cosparnh 

DATACOM 

FILENAME 

RANDCOM 

TERHGT 

THERM 

WAVEDAT 


MarsGRAM, 

MarsGRAM, 

MarsGRAM, 

MarsGRAM, 

MarsGRAM, 

MarsGRAM, 

MarsGRAM, 


cospar, Datastep, SETUP 

Block Data, Dustfact, Datastep, 

Datastep, SETUP 

Random, SETUP 

Terrain, SETUP 

ATMOS 2 , STEWART2 , SETUP 

Block Data, Wavepert 


SETUP 
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Table 2-3. Variables in the Common Blocks (Page 1 of 2) 


COSPARNH 



zc (159) 

REAL *4 

COSPAR heights (km) 

tc (159) 

REALM 

COSPAR temperatures (K) 

pc ( 159) 

REALM 

COSPAR pressures 

dc ( 1 5 9 ) 

REALM 

COSPAR densities 


DATACOM 

DTR 

REAL* 4 

factor for degrees to radians 

BETA 

REAL* 4 

parameter for computing viscosity of C02 

SVAL 

REAL* 4 

parameter for computing viscosity of C02 

DAY 

REAL* 4 

length of sidereal day (hours) 

CORFAC 

REAL* 4 

Coriolis factor for computing winds 

NPOS 

INTEGER* 4 

maximum number of positions 

NVARX 

INTEGER* 4 

x-code for plotable output (see Table 2-3) 

NVARY 

INTEGER* 4 

y-code for plotable output (see Table 2-3) 

logscale 

INTEGER* 4 

units for pressure and density output 

dustlat 

REAL* 4 

(0 = MKS, 1= log-base-10, 2-% deviation 
from COSPAR) 

latitude of local scale dust storm 

dust Ion 

REAL* 4 

longitude of local scale dust storm 

dusthgt 

REALM 

height of local scale dust storm 

radmax 

REAL *4 

maximum radius of local scale dust storm 

Rref 

REAL* 4 

Mars radius (km) of Cain 6 . 1 mb reference 

modpert 

INTEGER* 4 

ellipse 

perturbation model number 

aisO 

REAL *4 

initial areocentric longitude of sun 

intens 

REAL *4 

dust storm intensity value 

iuO 

INTEGER* 4 

unit number for screen output 

iup 

INTEGER* 4 

unit number for LIST output (set=Q for no 

maxf iles 

INTEGER* 4 

LIST) 

maximum number of files that compiler can open 


FILENAME 




lstfl 

CHAR* 12 

name of 

LIST file 

outf 1 

CHAR* 12 

name of 

OUTPUT file 


RANDCOM 







IX 

INTEGER* 4 

intermediate 

variable 

for 

random 

numbers 

IY 

INTEGER* 4 

intermediate 

variable 

for 

random 

numbers 

IZ 

INTEGER* 4 

intermediate 

variable 

for 

random 

numbers 


TERHGT 


terrain height versus latitude and longitude 

th ( 1 9 , 1 9 ) 

REAL *4 

THERM 


10.7 cm solar flux at Earth position (1AU) 

F107 

REAL* 4 

stdl 

REAL* 4 

number of std. dev. for thermospheric 



variation 
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Table 2.3. Variables in the Common Blocks (Page 2 of 2) 

WAVEDAT 


ampclOO ( 12 ) 

REAL *4 

diurnal wave amplitude vs height, clear, 
lat=00 

phclOO (12) 

REAL *4 

diurnal wave phase vs height, clear, lat=00 

ampdlOO ( 12 ) 

REAL * 4 

diurnal wave amplitude vs height, dusty, 
lat=00 

phdlOO { 12 ) 

REAL *4 

diurnal wave phase vs height, dusty, lat=00 

ampd200 { 12 ) 

REAL *4 

semidiurnal wave amplitude vs height, clear, 
lat=00 

phd200 { 12 ) 

REAL *4 

semidiurnal wave phase vs height, clear, 
lat=00 

ampcl2 0 ( 12 ) 

REAL *4 

diurnal wave amplitude vs height, clear, 
lat=2 0 

phcl20 (12) 

REAL *4 

diurnal wave phase vs height, clear, lat=20 

ampdl2 0 (12) 

REAL *4 

diurnal wave amplitude vs height, dusty, 
lat=2 0 

phdl20 {12) 

REALM 

diurnal wave phase vs height, dusty, lat=20 

ampd2 2 0(12) 

REAL * 4 

semidiurnal wave amplitude vs height, clear, 
lat=20 

phd220 { 12 ) 

REAL *4 

semidiurnal wave phase vs height, clear, 
lat=20 

ampcl45 (12) 

REAL *4 

diurnal wave amplitude vs height, clear, 
lat=45 

phcl45 (12) 

REAL * 4 

diurnal wave phase vs height, clear, lat=45 

ampdl45 (12) 

REAL* 4 

diurnal wave amplitude vs height, dusty, 
lat=45 

phdl45 (12) 

REAL* 4 

diurnal wave phase vs height, dusty, lat=45 

ampd245 (12) 

REAL *4 

semidiurnal wave amplitude vs height, clear, 
lat=45 

phd245 (12) 

REAL *4 

semidiurnal wave phase vs height, clear, 
lat=45 
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3. NEW MARS-GRAM FEATURES 


3.1 Stability-Limited Mountain Wave Perturbations 

A mountain-wave density perturbation model has been part of Mars-GRAM from the 
beginning. Recently noted is that effects of atmospheric stability should limit amount of 
growth with height attained from the mountain wave model. 

Deviation of the convective stability limit on density perturbations incorporated into 
Mars-GRAM 3.34 is as follows: 

Let the total atmospheric density be p = p 0 + p', where p 0 is the mean value and p' is 
the mountain-wave perturbed value. Static instability (density overturning) would result if the 
density were to increase with height. Thus stability requires that 


dp/dz = dpo/dz + dp'/dz < 0 (1) 

or 

dp'/dz < -dpo/dz = po/H , (2) 


where H is the scale height. If we approximate (dp'/dz) as 2 p' max /L, where L is the vertical 
scale (or wavelength), then the stability limit for density overturning is 


P'max/pO < L/(2H) 


(3) 


A more stringent stability condition (usually) than density overturning is the convective 
stability limit that temperature lapse rate (-dT/dz) not exceed adiabatic lapse rate (g/Cp, where 
g = gravity, Cp = specific heat at constant pressure). Thus, the convective stability constraint 
means that 


-(dTo/dz + dT’/dz) < g/Cp 


(4) 


or 


-dT’/dz < (g/Cp + dTo/dz) 


(5) 



If we approximate -dT'/dz as 2 T'^ / L , the convective stability limitation then becomes 


T'max / T 0 < (g/Cp + dTo/dz)L/(2T 0 ) , (6) 


or, in terms of the Brunt-Vaisala frequency, ©b, 


T'max / T 0 < COfi 2 L / (2 g) 


(7) 


If first order versions of the perturbed perfect gas law and hydrostatic balance relations, 


p7p 0 = p'/po + T’/To 

and 

d(p7p<>)/dz = (1/H)(T7T 0 ) 


( 8 ) 

(9) 


are used to convert the temperature perturbation in equation (7) to a density perturbation. The 
convective stability limit on density perturbations then becomes 


p7p 0 < [ cob 2 L / (2 g) ] [ 1 +L/(2H)]. (10) 


This convective stability limit, via equation (10), on allowable magnitude for mountain-wave 
perturbations has now been incorporated into the subroutine DATASTEP (DSTP199c), with 
a constraint that the density overturning condition of equation (3) also be met (DSTP199d). 
See further discussion of the DATASTEP subroutine in Appendix A. 


3.2 Comparison to COSPAR Reference Atmosphere 

Mars-GRAM Version 3.34 includes a comparison of model values of density, 
temperature, and pressure with those of the COSPAR reference atmosphere (Table XI of 
Pitts, et al., 1990?). COSPAR data values are read in as a height array, from a file named 
COSPAR.DAT. COSPAR values for a desired height are found by interpolating between 
heights in the tabular data. See discussion of the subroutine COSPAR in Appendix A for 
technical details of the methodology. 
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3.3 New Surface Temperature Parameterizations 

In early versions of Mars-GRAM, surface temperatures were estimated (TSURFACE 
subroutine) by empirically-derived regressions of temperature versus A, the daily average 
amount of solar radiation absorbed at the surface. Separate regressions were used for daily 
minimum and daily maximum temperature, via 


T = 

1 min 


hmin A + C m j n A 


and 


T m ax ~ <*max 




+ C n 


(ID 


( 12 ) 


with average daily temperature of (T min + T max )/2. 

A recent re-analysis of original Viking lander and InfraRed Thermal Mapper (IRTM) 
data, yielded an improved approach for estimating surface temperatures. First, an improved 
method is used to compute surface absorption, A, via 


A = < T > (1 - a)F 0 < cos(0) > 


(13) 


where < T > is the daily average solar transmittance, a is the surface albedo (from ALB 
function), F 0 is the top-of-atmosphere, direct-normal solar flux (for given latitude and day), 
and <cos(0)> is the daily average of cosine of solar zenith angle, 0. The average 
transmittance is computed by methods of Justus and Paris (1985), via 


<t> = 03/2 + ( 1 - GJ/2)exp(-S/po) , (14) 


where GJ is the single-scatter albedo of the dust (taken to be 0.85), 8 is the optical depth of the 
dust (taken to be 0.3), and po is the cosine of the noontime solar zenith angle, given by 


Po = sin( cp ) sin( <p s ) + cos( cp ) cos( (p s ) , (15) 

where <p is local latitude and (p s is latitude of the Sun. 

The simplified relation of equation (14) was compared with the results of a subroutine 
(FFACT) developed by Davies (1979) from accurate Monte Carlo radiative transfer 
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calculations. Equation (14) results were found to agree with FFACT values within a root- 
mean-square value of about 0.03, or roughly the same accuracy level as FFACT values 
reproduce original Monte Carlo simulations (Davies, 1979). 

With surface absorption, A, determined from equation (13), new regression relations 
were derived for T avg , the daily average temperature via 

T avg = a + bA + c A 2 , (16) 

if the latitude is outside the polar cap boundary, and 

T avg — ^cap b ca p A - c ca p P , (17) 

if latitude is inside the polar cap boundary. P is a polar cap correction that varies from 0 at 
the polar cap boundary to a maximum of 1 at the pole, when the cap boundary is at its largest 
seasonal value (factor polar in TSURFACE, line TSRF 35, Appendix A). 

New methodology for surface temperature estimates assumes that the daily range of 
surface temperatures, T amp = ( T max - T min ), is proportional to the daily range of surface 
absorption, Q, where Q is given by 

Q = <X> ( 1 - a ) F 0 [ < cos(0) > - H„ ] / 2 , (18) 

where the symbolism is the same as for equation (13) and p n is the midnight solar zenith 
angle, via 

|i n = sin( tp ) sin( (p s ) - cos( (p ) cos( <p s ) , (19) 

for latitudes at which the Sun is above the horizon at midnight; otherwise, |i n = 0. 

The daily surface temperature range is determined from 

T amp = 0.16 Q (20) 

and daily minimum and maximum surface temperatures are given by 
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( 21 ) 


T = T - T 

L nun 1 avg 1 amp 


Tmax — T avg + T amp 


( 22 ) 


For additional technical description of the new surface temperature methodology, see 
subroutine TSURFACE in Appendix A. 

The new regressions make relatively little difference in the average surface 
temperature. However, for cases in which the daily average absorption, A, is relatively large 
while the daily range in surface absorption, Q, is small, a significant change (reduction in 
daily temperature range) from the previous regressions can occur. This effect, illustrated in 
Figure 1 is most significant at high northern latitudes during northern hemisphere summer 
(near L s = 90°) and at high southern latitudes during southern hemisphere summer (near L s = 
270°). Changes in these seasons and latitude ranges would be apparent in revised plots of 
surface temperature analogous to Figures 1, 2, and 3 of Appendix B. The reduced diurnal 
range in temperatures that results from the new regressions is considered to be significantly 
more realistic in these cases. 
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Figure 3-1 . Comparison of Old and New Parameterizations for Surface Temperature in Mars- 
GRAM. The season is southern hemisphere summer (L s =270°). Lines are the new regression values 
for minimum, average and maximum daily surface temperature. Symbols are for the old regression 

results. 
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4. HOW TO RUN MARS-GRAM 


4.1 Program Input 

Two forms of Mars-GRAM are: interactive, in which values for all input options 
are provided interactively by the user at run time, and batch, in which values for all input 
options are provided by a NAMELIST format input file. Table 4-1 illustrates a sample 
operation of the interactive form of Mars-GRAM version 3.34 and Table 4-2 gives a 
sample of the NAMELIST file INPUT for the batch form of Mars-GRAM 3.34. 

For both the interactive or batch forms, values of the following input variables 
must be supplied: 


LSTFL 

Name of LIST file (see Table 4-4). For a listing to the console in 
interactive form enter filename CON. 

OUTFL 

Name of OUTPUT file (see Table 4-5). 

MONTH 

Month (1-12) for initial time 

MDAY 

Day of the month for initial time 

MYEAR 

Year for starting time is a 4-digit number. Alternative: input years 
1970-2069 as a 2-digit number 

NPOS 

Maximum number of positions to evaluate to automatically generate a 
profile. Use 0 if trajectory positions are read in from a TRAJDATA file. 

IHR 

Initial time, hour of the day GMT 

IMIN 

Initial time, minute of the hour 

SEC 

Initial time, second of the minute 

ALSO 

Value of the areocentric longitude of the Sun (Ls, in degrees) where a 
dust storm is to start. Use 0 if no dust storm is to be simulated. (Dust 
storm can be simulated only during the season of the Mars year for Ls 
between 180 and 320 degrees.) 

INTENS 

Dust storm intensity, an arbitrary intensity scale, with allowable values 
from 0.0 (no dust storm) to 3.0 (maximum intensity dust storm). [No 
prompt occurs in interactive form, if no dust storm is to be simulated 
(ALS0=0).] 



RADMAX 

DUSTLAT 

DUSTLON 

FI 07 

STDL 

MODPERT 

NR1 

NVARX 

NVARY 

LOGSCALE 


Maximum radius (km) a dust storm attains, developing according to the 
parameterized space and time profile of build up and decay in the 
program. If 0 or >10000 km is used, the storm is considered of global 
dimensions (uniformly covering the planet), but assumed to build up and 
decay in intensity according to the same temporal profile. [No prompt 
occurs in interactive, if no dust storm is to be simulated.] 

Latitude (degrees, North positive) for center of dust storm. [No prompt 
occurs in interactive, if no dust storm is to be simulated or storm has 
global dimensions.] 

Longitude (degrees, West positive) for center of dust storm. [No prompt 
occurs in interactive, if no dust storm is to be simulated or storm has 
global dimensions.] 

10.7-cm solar flux (units 10' 22 W/cm 2 ) at the average Earth orbit position 
(1 AU). Program automatically converts solar flux to its value at orbit 
position of Mars. 

Standard deviation parameter for long-term variations in the Stewart 
model thermosphere. Normal value is 0; allowable range is -3.0 to +3.0. 

Model number for perturbations to be computed: 1 is for random 
(mountain wave) model , 2 is for Zurek (tidal) wave model, and 3 
means use combined perturbations from both models. 

Seed value (integer) for random number generator. Allowable range is 1 
to 29999. No prompt occurs in interactive if MODPERT = 2. To do 
Monte Carlo simulations with a variety of perturbations, use a different 
random number seed on each run. To repeat a given perturbation 
sequence on a later run, use same random number seed value. 

x-code for plotable output (x-y pairs for 1-D line graphs or x-y-z triplets 
for 2-D contour plots). See Table 4-3 for list of variables associated 
with x-code (e.g., if NVARX = 1, output is for plotting versus height 
above reference ellipsoid). 

y-code for 2-D contour plot output (x-y-z triplets). Use y-code 0 for 1-D 
line graph (x-y pair) plots. See Table 4-3 for list of y-code values and 
parameters. 

Parameter controls units of output values for density and pressure on 
output plot files. Value 0 means use regular density and pressure units 
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(kg/m^ and N/m^); 1 means output logarithm (base 10) in regular units; 
and 2 means output percent deviation from COSPAR values. 

FLAT Latitude of initial point to simulate (degrees. North positive) 

FLON Longitude of initial point to simulate (degrees, West positive) 

FHGT Height (km) of initial point to simulate, above reference ellipsoid 

DELHGT Height increment (km) between successive steps in automatically 
generated profile (positive upward) 

DELLAT Latitude increment (degrees, Northward positive) between successive 
steps in automatically generated profile 

DELLON Longitude increment (degrees. Westward positive) between successive 
steps in automatically generated profile 

DELTIME Time increment (seconds) between steps in automatically generated 
profile 

Two auxiliary input files are also required. File HEIGHTS.DAT contains terrain 
height data array (terrain height, km, above the reference ellipsoid - see explanation in 
description of subroutine TERRAIN in Appendix A). File COSPAR.DAT contains height 
profile of COSPAR temperature, density, and pressure values (see explanation in 
description of COSPAR subroutine in Appendix A). 

If the pre-computed trajectory mode is used (NPOS=0), read trajectory data from 
TRAJDATA file. Each line of TRAJDATA file is a position and time to compute 
atmospheric parameters. Input lines contain time (seconds, from initial time), height (km, 
relative to reference ellipsoid), latitude (degrees, North positive), and longitude (degrees. 
West positive). 

For automatically-generated profiles, output continues until the maximum number 
of positions (NPOS) is reached. For trajectory positions, enter input from TRAJDATA 
file, output continues until end of the file is reached. For interactive, the program prompts 
for additional input values for initial date and number of positions. The program is 
terminated by giving values of 0 for requested input. (See end of Table 4-1). 
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Table 4-1. Sample Operation of the Interactive Form of Mars-GRAM 3.34 

Mars-GRAM Interactive version 3.34 - November 1, 1995 
Enter name for LIST file (CON for console listing) : 

VIKING1 . LST 

Enter name for OUTPUT file: 

VIKINGl . OUT 

Enter Month, Day of Month, 4-digit Year, and Max Number Positions 
7 20 76 21 

Enter initial GMT Time (Hours, Minutes, Seconds) 

12 30 0 

Ls = 97.0 degrees for this date. 

Dust storms can occur between Ls = 180 and Ls = 320. 

Enter starting Ls value for dust storm (or 0 for none) . 

0 

Enter mean F10.7 flux at 1AU (nominal value = 150) 

and +/- number of std. deviations for thermosphere variation 

185 0 

Enter perturbation model: l=random, 2=wave, 3=both 

Enter Starting Random Number (any positive integer < 30,000) 

1001 . J 

x— code and y— code for plotable output versus desired parameter (s) . 


Code Parameter 


1 Height (above reference ellipsoid, km) 

2 Height (above local terrain, km) 

3 Latitude (deg.) 

4 West Longitude (deg.) 

5 Time from start (Earth seconds) 

6 Time from start (Martian Sols) 

7 Areocentric Longitude of Sun, Ls (deg.) 

8 Hour Angle for Local Time (Mars hours * 15) 

Use y-code = 0 for plotable output vs x-code variable only 
2 0 

For density and pressure data units, enter: 

0 for normal (MKS) , 1 for log-base-10, 2 for % difference from COSPAR 

2 

Enter Initial Latitude (deg.), West Longitude (deg.) 

22 48 

Surface elevation = -.50 km at this location 

Enter Initial Height (km) 

-0.5 

Enter Increments in Height (km), Latitude (deg.), West Longitude (deg.), 
and Time (sec.) 

10 0 0 0 
Computing data. 

Enter Month, Day of Month, 4-digit Year, and Max Number Positions 
0 0 0 0 

Normal Termination 
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Table 4-2. Sample NAMELIST file INPUT for Batch Form of Mars-GRAM 3.34. 
Inline comments are appended after the ! symbol. Some FORTRAN compilers do 
not allow inline comments in NAMELIST data files, in which case they cannot be 
included in the INPUT. Some FORTRAN compilers use alternate forms of the initial 
and final lines of the file (e.g. &INPUT for $INPUT and / for $END is the form for 
Microsoft DOS FORTRAN). 


$ INPUT 

LSTFL = ' VIKINGl . LST ' , ! List file name (CON for console listing) 

OUTFL = * VIKINGl . OUT ' , ! Output file name 

MONTH = 7, ! month of year 

MDAY =20, ! day of month 

MYEAR = 76, ! year (4-digit; 1970-2069 can be 2-digit) 

NPOS =21, ! max # positions to evaluate (0 = read data from file) 

IHR =12, ! GMT hour of day 

IMIN =30, ! minute of hour 

SEC =0.0, ! second of minute (for initial position) 

ALSO =0.0, ! starting Ls value (degrees) for dust storm (0 = none) 

INTENS =0.0, ! dust storm intensity (0.0 - 3.0) 

RADMAX =0.0, ! max. radius (km) of dust storm (0 or >10000 = global) 

DUSTLAT =0.0, ! latitude (deg) for center of dust storm 

DUSTLON =0.0, ! West longitude (deg) for center of dust storm 

F107 = 185.0, ! 10.7 cm solar flux (10**-22 W/cm**2, at 1 AU) 

STDL =0.0, ! std. dev. for thermosphere variation (-3.0 to +3.0) 

MODPERT =3, ! perturbation model; l=random, 2=wave, 3=both 

NRl =1001, ! starting random number (0 < NR1 < 30000) 

NVARX =2, ! x-code for plotable output (l=hgt above ref. ellipse) 

NVARY =0, ! y-code for 2-D plotable output (0 for 1-D plots) 

LOGSCALE = 0, ! 0=regular density, l=log (density) , 2=COSPAR deviations 

FLAT = 22.0, ! initial latitude (N positive), degrees 

FLON =48.0, ! initial longitude (West positive), degrees 

FHGT = -0.5, ! initial height (km), above ref. ellipse 

DELHGT =10.0, ! height increment (km) between steps 

DELLAT =0.0, ! latitude increment (deg) between steps 

DELLON =0.0, ! West longitude increment (deg) between steps 

DELTIME =0.0, ! time increment (sec) between steps 


$END 
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Table 4-3. List of x-code and y-code Values and Parameters for Plotable Output 

(graphics) Files 


Code Parameter SETU283 

SETU284 

1 Height (above reference ellipsoid, km) SETU285 

2 Height (above local terrain, km) SETU286 

3 Latitude (deg.) SETU287 

4 West Longitude (deg.) SETU288 

5 Time from start (Earth seconds) SETU289 

6 Time from start (Martian Sols) SETU290 

7 Areocentric Longitude of Sun, Ls (deg.) SETU291 

8 Hour Angle for Local Time (Mars hours * 15) SETU292 


Use y-code = 0 for plotable output vs x-code variable only SETU293 
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4.2 Program Output 

Three general types of program output are: (1) a "LIST" file, containing header and 
descriptor information, suitable for printing or viewing by an analyst (LIST file, Table 4- 
4), (2) "OUTPUT" file, containing no header or descriptor information, one line per 
output position, suitable for reading into another program for additional analysis 
(OUTPUT file, in Table 4-5), and (3) set of "plotable" output files, suitable for input to a 
graphics program (Table 4-6). 

Plot output files contain either x-y data pairs or x-y-z data triplets, determined from 
selected values for x-code (NVARX) and y-code (NVARY). If 1-D line-graph (x-y pair) 
data is the selected plot output option (y-code = 0), then the x-code variable normally 
appears first in the output pair. However, if the x-code variable selected for plot output is 
a height variable (x-code = 1 or 2) then the plot output values give the height as the second 
variable of the output pair. This facilitates graphs with height as the ordinate, since many 
graphics programs expect input data pairs in abscissa-ordinate order. If 2-D contour plot 
(x-y-z triplet) data is the selected plot output option (y-code * 0), then the plot output 
variables always appear as: x-code variable, y-code variable, plot output variable. See 
Table 4-7 for list of plot output variables generated and associated file unit numbers (plot 
output files are unit numbers 2 1 through 28 and 30 through 34). 

Table 4-7 gives the unit numbers for all program input and output files. With the 
exception of the screen output and the TRAJDATA input file, the program input and 
output files are also referred to within the program by a FILES array. The FILES array 
index values associated with the program input and output files are also shown in Table 4- 
7. 


If the user desires to suppress the LIST, OUTPUT, and plotable output files (to 
handle output in a user-provided program), set the LIST file unit number (iup) to 0 (line 
BLKD 15, in BLOCK DATA routine in MARSSUBS.FOR file) and re-compile the 
program. The unit number associated with "screen" output (iuO), normally 6, is set to any 
other value by changing it at line BLKD 16 and re -compiling the program. 

WRITE statements and output FORMATS used to produce the OUTPUT file are 
shown in the following code lines from the DATASTEP subroutine (in MARSSUBS.FOR 
file): 


If (logscale.eq.OjThen DSTP34 6b 

WRI TE ( 2 9 , 8 0 0 ) C SEC , VAR , CLAT , CLON , dens , TEMP , EWWIND , DSTP 3 4 7 

Sl NSWIND, SIGD, DENSWA,ALS DSTP348 

Else DSTP348& 

WRITE (29,810) CSEC , VAR, CLAT, CLON, dens , TEMP, EWWIND, DSTP 3 4 8b 

& NSWIND, SIGD, DENSWA, ALS DSTP348c 

Endif DSTP348d 

Endif DSTP348e 
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800 FORMAT (F10 .0, 3F7 . 2 , 1PE9 . 2 , 0F3F7 .1, 3F6 . 1) 
8X0 FORMAT ( F10 .0, 3F7 . 2 , F9 . 2 , 3F7 . 1, 3F6 . 1 ) 


DSTP349 

DSTF349* 


If other variables or output formats are desired, alter appropriate program lines and re- 
compile the program. 

The output data, some echoed from input values for a record of the options 
selected, are illustrated by the sample LIST file in Table 4-4. Output variables and 
descriptions are as follows: 

MONTH, MDAY, MYEAR, IHR, IMIN, SEC, ALSO, INTENS, RADMAX, DUSTLAT, 
DUSTLON, FI 07, STDL, MODPERT, and NR1 are defined in Section 4.1. 


DATE 

CSEC 

ALS 

OHGT 

OHGTS 

HSCALE 

CLAT 

CLON 

SUNLON 

TLOCAL 

TEMP 

PRES 

DENSLO 


Julian date, computed from input year, month and day 

Time (seconds) from start of simulation [also units of Mars days (sols)] 

Areocentric longitude of the Sun (Ls) in degrees (angle at Mars analogous to 
right ascension of the Sun for the Earth) 

Height above reference ellipsoid (km) 

Height above local terrain (km) (shown in parentheses in LIST file output) 
Pressure scale height (km) 

Current latitude (positive North) 

Current longitude (positive West) 

Mars longitude of sub-solar point (degrees) 

Local solar time (in Mars hours, 1 /24th sols) 

Temperature (K) at current location (also output to plot file TEMP, unit 25, 
Table 4-7) 

Pressure (N/m^) at current location (also output to plot file PRES, unit 26, 
Table 4-7; units in plot file controlled by LOGSCALE option) 

Low (approximately -1 standard deviation) density (kg/m^) at current 
location (also output to plot file DENSLO, unit 33, Table 4-7 ; units in plot 
file controlled by LOGSCALE option) 



DENS Average density (kg/m^) at current location (also output to plot file 

DENS A V, unit 22, Table 4-7; units in plot file controlled by LOGSCALE 
option) 

DENSHI High (approximately +1 standard deviation) density (kg/m 3) at current 

location (also output to plot file DENSHI, unit 34, Table 4-7; units in plot 
file controlled by LOGSCALE option) 

DEVLO Low density expressed as percentage deviation from COSPAR value 

DEVAV Average density expressed as percentage deviation from COSPAR value 

DEVHI High density expressed as percentage deviation from COSPAR value 

DENSP Density perturbation (relative to mean value), expressed as percentage of 
mean (also output to plot file DENSRP, unit 24, Table 4-7). 

EWWIND Eastward wind component (m/s) (also output to plot file EWWIND, unit 27, 
Table 4-7). 

NSWIND Northward wind component (m/s) (also output to plot file NSWIND, unit 28, 
Table 4-7). 

Other variables not output to LIST file but output to plot files (may be selected for 

output to OUTPUT file) are the following: 

SIGD Standard deviation of density perturbations (percent of mean value) (output 
to plot file DENSRM, unit 21, Table 4-7) 

DENSWA Value of the Zurek wave model density perturbation (percent of mean value) 
(output to plot file DENSWA, unit 23, Table 4-7) 

TMAX Daily maximum temperature (K) at surface, for current location (output to 
plot file TMAX, unit 30, Table 4-7) 

TMIN Daily minimum temperature (K) at surface, for current location (output to 
plot file TMIN, unit 31, Table 4-7) 

TAVG Daily average temperature (K) at surface, for current location (output to plot 
file TAVG, unit 32, Table 4-7) 
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Table 4-4. LIST file (VIKING1.LST) Produced by Either Interactive or Batch Form 

(Page 1 of 4) 


Mars-GRAM Version 3.34 Interactive Form - November 1, 1995 

Date = 7/20/1976 Julian Date = 2442980.0 GMT Time * 12:30: .0 

F10.7 flux = 185.0 (1 AU) 68.1 (Mars), standard deviation = .0 

Perturbation model = 3 Starting random number = 1001 

Time (rel. to TO) * .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = -.50 km ( .00 km) Scale Height = 12.48 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 243.4 K Pressure = 7.478E+02 N/m**2 

Density (Low, Avg. , High) = 1.502E-02 1.607E-02 1.712E-02 kg/m**3 

Departure, COSPAR NH Mean = -7.5% -1.0% 5.5% 

Density perturbation = 8.76 % of mean value 

Eastward Wind = 4.3 m/s Northward Wind = -.3 m/s 


Time (rel . to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 9.50 km ( 10.00 km) Scale Height = 10.39 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 202.7 K Pressure = 3.055E+02 N/m**2 

Density {Low, Avg., High) = 7.619E-03 7.884E-03 8.152E-03 kg/m**3 

Departure, COSPAR NH Mean = 12.9 % 16.8 % 20.8 % 

Density perturbation = -1.30 % of mean value 

Eastward Wind = 1.3 m/s Northward Wind = -5.5 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 19.50 km ( 20.00 km) Scale Height = 9.49 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 185.1 K Pressure = 1.122E+02 N/m**2 

Density (Low, Avg., High) = 3.038E-03 3.172E-03 3.309E-03 kg/m**3 

Departure, COSPAR NH Mean = 10.2 % 15.0 % 20.0 % 

Density perturbation = -1.41 % of mean value 

Eastward Wind = -1.8 m/s Northward Wind = -9.0 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 29.50 km ( 30.00 km) Scale Height = 8.72 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 170.2 K Pressure = 3.800E+01 N/m**2 

Density (Low, Avg. , High) = 1.096E-03 1.168E-03 1.242E-03 kg/m**3 

Departure, COSPAR NH Mean = 6.0% 13.0% 20.1% 

Density perturbation = -.30 % of mean value 

Eastward Wind = -4.8 m/s Northward Wind = -13.9 m/ s 


Time (rel. to TO) = .0 sec. { .000 sols) Ls = 97.0 deg. 

Height = 39.50 km { 40.00 km) Scale Height = 8.19 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 159.8 K Pressure = 1.192E+01 N/m**2 

Density (Low, Avg., High) = 3.655E-04 3.903E-04 4.160E-04 kg/m**3 

Departure, COSPAR NH Mean = 1.8% 8.7% 15.9% 

Density perturbation = -.32 % of mean value 

Eastward Wind = -8.2 m/s Northward Wind = -20.1 m/s 


Time (rel. to TO) = -0 sec. ( .000 sols) Ls - 97.0 deg. 

Height = 49.50 km { 50.00 km) Scale Height = 7 . 67 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 149.6 K Pressure = 3.488E+00 N/m**2 

Density (Low, Avg., High) = 1.103E-04 1.220E-04 1.340E-04 kg/m**3 

Departure, COSPAR NH Mean = -4.0% 6.3% 16.8% 
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Table 4-4. LIST file (VIKING1.LST) Produced by Either Interactive or Batch Form 

(Page 2 of 4) 

Density perturbation = 5.36 % of mean value 

Eastward Wind = -11.9 m/s Northward Wind = -26.7 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 59.50 km { 60.00 km) Scale Height = 7.39 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 144.1 K Pressure = 9.618E-01 N/m**2 

Density (Low, Avg., High) = 3.018E-05 3.491E-05 3.973E-05 kg/m**3 

Departure, COSPAR NH Mean = -10.9% 3.0% 17.3% 

Density perturbation = -9.27 % of mean value 

Eastward Wind = -16.2 m/s Northward Wind = -33.9 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 69.50 km ( 70.00 km) Scale Height = 7.13 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 139.0 K Pressure = 2.552E-01 N/m**2 

Density (Low, Avg., High) = 8.298E-06 9.599E-06 1.092E-05 kg/m**3 

Departure, COSPAR NH Mean = -11.1 % 2.9 % 17 . 1 % 

Density perturbation = -1.71 % of mean value 

Eastward Wind = -20.7 m/s Northward Wind = -39.8 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 79.50 km ( 80.00 km) Scale Height = 7.34 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 136.5 K Pressure = 5.487E-02 N/m**2 

Density (Low, Avg., High) = 1.680E-06 2.100E-06 2.549E-06 kg/m**3 

Departure, COSPAR NH Mean = -31.4 % -14.3 % 4.1 % 

Density perturbation = 3.65 % of mean value 

Eastward Wind = -23.7 m/s Northward Wind = -40.5 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 89.50 km ( 90.00 km) Scale Height = 7.55 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 139.2 K Pressure = 8.699E-03 N/m**2 

Density (Low, Avg., High) = 2.227E-07 3.256E-07 4.505E-07 kg/m**3 

Departure, COSPAR NH Mean = -65.3 % -49.3 % -29.9 % 

Density perturbation = -31.47 % of mean value 

Eastward Wind = -22.0 m/s Northward Wind = -31.5 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 99.50 km ( 100.00 km) Scale Height = 7.90 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 144.4 K Pressure = 1.470E-03 N/m**2 

Density (Low, Avg., High) = 3.083E-08 5.284E-08 8.334E-08 kg/m**3 

Departure, COSPAR NH Mean = -81.9 % -68.9 % -51.0 % 

Density perturbation = -7.08 % of mean value 

Eastward Wind = -18.3 m/s Northward Wind = -21.2 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 109.50 km ( 110.00 km) Scale Height = 8.44 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 152.9 K Pressure = 2.776E-04 N/m**2 

Density (Low, Avg., High) = 4.922E-09 9.394E-09 1.618E-08 kg/m**3 

Departure, COSPAR NH Mean = -88.9 % -78.7 % -63.4 % 

Density perturbation = 128.57 % of mean value 

Eastward Wind = -11.5 m/s Northward Wind = -9.8 m/s 
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Table 4-4. LIST file (VIKING1.LST) Produced by Either Interactive or Batch Form 

(Page 3 of 4) 

Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 119.50 km ( 120.00 km) Scale Height = 11.18 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 199.8 K Pressure = 7.618E-05 N/m**2 

Density (Low, Avg., High) = 1.317E-09 1.959E-09 2.886E-09 kg/m**3 

Departure, COSPAR NH Mean = -89.6 % -84.5 % -77.2 % 

Density perturbation = -17.49 % of mean value 

Eastward Wind = -4.4 m/s Northward Wind = -3.2 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 129.50 km ( 130.00 km) Scale Height = 14.53 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 255.1 K Pressure = 3.512E-05 N/m**2 

Density (Low, Avg., High) = 4.979E-10 6.990E-10 9.288E-10 kg/m**3 

Departure, COSPAR NH Mean = -87.5 % -82.4 % -76.6 % 

Density perturbation = 1.15 % of mean value 

Eastward Wind = -2.0 m/s Northward Wind = -2.8 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 139.50 km ( 140.00 km) Scale Height = 16.60 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 285.8 K Pressure = 1.851E-05 N/m**2 

Density (Low, Avg., High) = 2.262E-10 3.241E-10 4.298E-10 kg/m**3 

Departure, COSPAR NH Mean = -80.3 % -71.8 % -62.6 % 

Density perturbation = 1.73 % of mean value 

Eastward Wind = -.8 m/s Northward Wind = -1.7 m/s 


Time (rel. to TO) = -0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 149.50 km ( 150.00 km) Scale Height = 18.01 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 302.9 K Pressure = 1.039E-05 N/m**2 

Density (Low, Avg., High) = 1.124E-10 1.687E-10 2.300E-10 kg/m**3 

Departure, COSPAR NH Mean = -77.1 % -65.7 % -53.2 % 

Density perturbation = 37.30 % of mean value 

Eastward Wind = -.4 m/ s Northward Wind = -1.0 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 159.50 km ( 160.00 km) Scale Height = 19.10 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 312.5 K Pressure = 6.064E-06 N/m**2 

Density (Low, Avg., High) = 5.869E-11 9.339E-11 1.324E-10 kg/m**3 

Departure, COSPAR NH Mean = -76.6 % -62.8 % -47.2 % 

Density perturbation = -35.75 % of mean value 

Eastward Wind = -.3 m/s Northward Wind = -.5 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 169.50 km ( 170.00 km) Scale Height = 20.07 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 317.9 K Pressure = 3.638E-06 N/m**2 

Density (Low, Avg., High) = 3.170E-11 5.362E-11 7.946E-11 kg/m**3 

Departure, COSPAR NH Mean = -77.2 % -61.4 % -42.8 % 

Density perturbation = -8.43 % of mean value 

Eastward Wind = -.2 m/s Northward Wind = -.3 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 179.50 km ( 180.00 km) Scale Height = 21.06 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 
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Table 4-4. LIST file ( VIKING l.LST) Produced by Either Interactive or Batch Form 

(Page 4 of 4) 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 321.0 K Pressure = 2.236E-06 N/m**2 

Density (Low, Avg., High) = 1.757E-11 3.158E-11 4.897E-11 kg/m**3 

Departure, COSPAR NH Mean = -78.3 % -61.1 % -39.6 % 

Density perturbation = -.26 % of mean value 

Eastward Wind = -.2 m/s Northward Wind = -.1 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 189.50 km ( 190.00 km) Scale Height = 22.15 km 
Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 322.7 K Pressure = 1.407E-06 N/m**2 

Density (Low, Avg., High) = 9.968E-12 1.899E-11 3.078E-11 kg/m**3 

Departure, COSPAR NH Mean = -79.5 % -60.9 % -36.7 % 

Density perturbation = -3.07 % of mean value 

Eastward Wind = -.1 m/s Northward Wind = -.1 m/s 


Time (rel. to TO) = .0 sec. ( .000 sols) Ls = 97.0 deg. 

Height = 199.50 km ( 200.00 km) Scale Height = 23.43 km 

Latitude = 22.000 degrees West Longitude = 48.000 degrees 

Sun Longitude = 111.188 deg. Local Time = 16.21 Mars hours 

Temperature = 323.7 K Pressure = 9.066E-07 N/m**2 

Density (Low, Avg., High) = 5.795E-12 1.164E-11 1.967E-11 kg/m**3 

Departure, COSPAR NH Mean = -80.6 % -61.1 % -34.3 % 

Density perturbation = -4.85 % of mean value 

Eastward Wind = -.1 m/s Northward Wind = .0 m/s 
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Table 4-5. OUTPUT file (VIKINGl.OUT) Produced by Either Interactive or Batch 

Form (Heading added for readability) 


Dens 


Time 

Hgt 

Lat 

Lon 

%COSPAR 

0 . 


00 

22 . 

00 

48 . 

00 

- 1 . 

00 

0 . 

10 . 

00 

22 . 

00 

48 . 

00 

16 . 

80 

0 . 

20 . 

00 

22 . 

00 

48 . 

00 

15 . 

05 

0 . 

30 . 

00 

22 . 

00 

48 . 

00 

13 . 

01 

0 . 

40 . 

00 

22 . 

00 

48 . 

00 

8 . 

71 

0 . 

50 . 

00 

22 . 

,00 

48 . 

00 

6 . 

27 

0 . 

60 . 

00 

22 . 

, 00 

48 . 

.00 

3 . 

04 

0 . 

70 . 

.00 

22 . 

.00 

48 . 

,00 

2 . 

89 

0 . 

80 . 

00 

22 . 

.00 

48 . 

.00 

- 14 . 

.28 

0 . 

90 . 

.00 

22 . 

.00 

48 . 

.00 

- 49 . 

.33 

0 . 

100 . 

.00 

22 , 

.00 

48 . 

.00 

- 68 . 

.94 

0 . 

110 , 

.00 

22 

.00 

48 . 

.00 

- 78 . 

.75 

0 . 

120 . 

.00 

22 

.00 

48 . 

.00 

- 84 . 

. 54 

0 . 

130 . 

.00 

22 

.00 

48 

.00 

-82 , 

.41 

0 . 

140 , 

.00 

22 

.00 

48 

.00 

-71 

.78 

0 . 

150 

.00 

22 

.00 

48 

.00 

-65 

.68 

0 . 

160 

.00 

22 

.00 

48 

. 00 

-62 

.78 

0 . 

170 

.00 

22 

.00 

48 

.00 

-61 

.41 

0 . 

180 

.00 

22 

.00 

48 

.00 

-61 

.06 

0 . 

190 

.00 

22 

.00 

48 

.00 

-60 

.94 

0 . 

200 

.00 

22 

.00 

48 

.00 

-61 

.10 


Temp 

EWind 

Nwind 

Dev 

Wave 

Ls 

243 . 

4 

4 . 

3 


3 


9 

5 . 

6 

97 . 

0 

202 . 

7 

1 . 

3 

- 5 . 

5 

1 . 

9 

1 . 

5 

97 . 

0 

185 . 

1 

- 1 . 

,8 

- 9 . 

0 

2 . 

9 

1 . 

4 

97 . 

0 

170 . 

2 

- 4 . 

.8 

-13 . 

9 

3 . 

9 

2 . 

4 

97 . 

0 

159 . 

,8 

- 8 . 

.2 

- 20 . 

, 1 

4 . 

8 

1 . 

6 

97 . 

0 

149 . 

,6 

- 11 . 

,9 

- 26 . 

.7 

4 . 

9 

4 . 

9 

97 . 

0 

144 . 

.1 

- 16 . 

.2 

-33 . 

.9 

4 . 

.9 

8 . 

,8 

97 . 

0 

139 . 

.0 

- 20 . 

.7 

- 39 . 

.8 

4 . 

.9 

8 . 

,8 

97 . 

0 

136 , 

.5 

-23 . 

.7 

- 40 . 

.5 

11 . 

.9 

8 . 

, 8 

97 . 

0 

139 . 

.2 

- 22 . 

.0 

- 31 , 

.5 

26 . 

.2 

8 . 

.8 

97 . 

0 

144 . 

.4 

- 18 . 

.3 

-21 

.2 

40 . 

.9 

8 . 

, 8 

97 . 

0 

152 , 

.9 

- 11 . 

.5 

-9 

.8 

51 . 

. 1 

8 . 

. 8 

97 . 

,0 

199 , 

.8 

-4 

.4 

-3 

.2 

31 

.2 

8 . 

.8 

97 . 

,0 

255 

.1 

-2 

.0 

-2 

.8 

22 

.0 

8 . 

.8 

97 . 

.0 

285 

.8 

- 

.8 

-1 

.7 

22 

.6 

8 

.8 

97 . 

.0 

302 

.9 

- 

.4 

-1 

.0 

26 

.1 

8 

.8 

97 . 

.0 

312 

.5 

- 

.3 

- 

.5 

30 

.7 

8 

.8 

97 . 

. 0 

317 

.9 

- 

.2 

- 

.3 

35 

.7 

8 

. 8 

97 . 

.0 

321 

.0 

- 

.2 

- 

.1 

40 

.9 

8 

.8 

97 . 

.0 

322 

.7 

- 

. 1 

- 

.1 

46 

.0 

8 

. 8 

97 

.0 

323 

.7 

- 

.1 


.0 

50 

. 8 

8 

.8 

97 

.0 
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Table 4-6. Sample Plotable Output File (DENSAV) [first variable = average density 
(% COSPAR difference), second variable = height] 


-.9976 

. 0000 

16 

. 80 

10 . 

00 

15 

. 05 

20 . 

00 

13 

.01 

30 . 

00 

8 . 

714 

40 . 

00 

6 . 

274 

50 . 

00 

3 . 

039 

60 . 

00 

2 . 

891 

70 . 

00 

-14 

.28 

80 . 

00 

-49 

.33 

90 . 

00 

-68 

.94 

100 

.0 

-78 

.75 

110 

.0 

-84 

.54 

120 

. 0 

-82 

.41 

130 

. 0 

-71 

.78 

140 

.0 

-65 

.68 

150 

.0 

-62 

.78 

160 

. 0 

-61 

.41 

170 

.0 

-61 

.06 

180 

.0 

-60 

. 94 

190 

.0 

-61 

. 10 

200 

.0 
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Table 4-7. File Names Used in the Mars-GRAM Programs 


FILES 

Array 

Index 


File Name (description of file) 


Unit Number 


1 

LSTFL 

(default name = LIST) 

iup 

Set 

and 

2 

DENSRM 

(std. deviation in density) 

21 

3 

DENSAV 

(average density) 

22 

4 

DENSWA 

(wave model perturbations) 

23 

5 

DENSRP 

(random perturbations) 

24 

6 

TEMP 

(temperatures) 

25 

7 

PRES 

(pressure) 

26 

8 

EWWIND 

(eastward wind component) 

27 

9 

NSWIND 

(northward wind component) 

28 

10 

OUTFL 

(default name = OUTPUT) 

29 

11 

HEIGHTS 

.DAT (terrain heights) 

9 

12 

TMAX 

(daily max. surf, temp.) 

30 

13 

TMIN 

(daily min. surf, temp.) 

31 

14 

TAVG 

(daily avg. surf, temp.) 

32 

15 

DENSLO 

(low density, -1 sigma) 

33 

16 

DENSHI 

(high density, +1 sigma) 

34 

17 

COSPAR. 

DAT (COSPAR data arrays) 

10 

N/A 

TRAJDATA (trajectory input) 

7 

N/A 

Unnamed 

l screen output 

iuO 


0 to suppress LIST 


(default=6 BLKD 16) 
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4.3 How to Use Mars-GRAM Batch Form as Subroutines in Other 
Programs 


Stripped of all embedded comment statements, the batch form program 
(MARSGRMB . FOR) consists of the following lines of code: 


DOUBLE PRECISION DATEO 
Real NSWIND 
Integer EOF 

Call Setup (CHGT, CLAT, CLON, CSEC, DATEO , RHO, DELHGT, DELLAT, DELLON, 
& DELTIME , MAXNUM ) 

DO 900 I ■ 0, MAXNUM 

Call Data8tep ( I , CHGT, CLAT, CLON, CSEC, DATEO , RHO, EOF, DELHGT, 

& DELLAT , DELLON , DELTIME , TEMP , PRES , DENSLO , DENS , DENSHI , DENSP , 

& EWWIND , NSWIND ) 

If (EOF .eg. l)Goto 999 
900 Continue 

999 STOP * Normal Termination' 

END 


All variables are initialized through SETUP subroutine which must be called once 
at beginning of program. Atmospheric parameters are evaluated at successive 
locations (either trajectory or automatic profile mode) with DATASTEP 
subroutine. Any desired output (or transfer to other subroutines) is done 
within the loop that calls DATASTEP (i.e. just after "Call Datastep" statement). 
To suppress the normal LIST, OUTPUT, and graphics file output, set iup=0 at line 
BLKD 15 in the BLOCK DATA routine (in the MARSSUBS.FOR file) . 

To embed SETUP and DATASTEP subroutines into a program that computes 
trajectory or orbit positions, the calling program keeps current position 
updated (height CHGT, latitude CLAT, longitude CLON, at time after start of the 
run CSEC) and computes velocity components from which displacements in height, 
latitude, and longitude can be computed for any desired time step. 

A program using SETUP and DATASTEP as subroutines takes the following form 
(added pseudo-code is shown in brackets) : 

DOUBLE PRECISION DATEO 
Real NSWIND 
Integer EOF 

Call Setup (CHGT, CLAT, CLON, CSEC, DATEO , RHO, DELHGT, DELLAT, DELLON, 

& DELTIME, MAXNUM) 

[Comment: Evaluate the atmospheric parameters at the initial 
position by calling Datastep with 1*0, so that the current 
position is not updated.] 

Call Datastep ( 0 , CHGT, CLAT, CLON, CSEC, DATEO, RHO, EOF, DELHGT, 

& DELLAT, DELLON, DELTIME, TEMP, PRES, DENSLO, DENS, DENSHI, DENSP, 

& EWWIND, NSWIND) 

100 Continue 

[Use the vertical velocity in the trajectory to compute the 
height displacement, DELHGT, the north-south velocity to 
compute the latitude displacement, DELLAT, and the east -west 
velocity to compute the longitude displacement, DELLON, based 
on the desired time displacement, DELTIME. These displacement 
values will override the values initially defined in the 
SETUP routine] 


Call Datastep (1, CHGT, CLAT, CLON, CSEC, DATEO, RHO, EOF, DELHGT, 
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& DELLAT , DELLON, DELTIME , TEMP , PRES , DENSLO , DENS , DENSHI , DENSP , 

& EWWIND , NSWIND ) 

[Comments With I * 0 in the call to Datastep, the position 
is updated by the position displacements, before the new 
atmospheric parameters are evaluated.] 

[Write out the current position (CHOT, CLAT, CLON) , as it 
has been updated by Datastep, along with any other output 
desired. Any other desired analysis of the variables at the 
current position can also be done here in the program. ] 

[If the last position and time has been computed, go to 999. 
Otherwise, go to 100.] 

999 STOP ' Normal Termination' 

END 
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5. DIAGNOSTIC AND PROGRESS MESSAGES 
5.1 Interactive Form Main Program Mars-GRAM 

The interactive form Mars-GRAM main program prompts the user for input 
options. Some prompt statements provide legal values that can be acceptable for the 
prompted input option. If the program repeats the prompt message after entered input, a 
value outside the legal range was input. Corrective action is to input a value within the 
legal range. Prompt message ceases and operation continues. 

Other specific diagnostic or progress messages follow. All STOP messages and 
messages written to unit iuO normally appear on the interactive user screen (default value 
for iuO is 6, set at line BLKD 16 in BLOCK DATA). 


Code Producing the Message: 

50 Writ© ( iuO , 60 ) files ( j ) , ierr ( j ) MARS108 

€0 Format ( lx, a 12 , ' File open error! Error =',i5) MARS109 

Purpose: 

To indicate an error occurred in attempt to open one of the 17 files in FILES array 
(see Table 4-7 for file names). Message gives file name and error code that caused the 
(first encountered) problem. 

Remedy: 

Depending on meaning of error code on compiler system, take appropriate 
corrective action to ensure proper file open operation. Rerun program. 


Code Producing the Message: 

65 Stop ' Error or EOF on COSPAR.DAT file!' 


MARS 1 12 g 


Purpose: 



To indicate an End-of-File or read error occurred on COSPAR.DAT file (file may 
not exist or path to this file may not be set appropriately). 

Remedy: 

Correct error condition (i.e., copy or restore file to right location). Rerun program. 


Code Producing the Message: 

80 If (lat .na. -100+10*1) Stop ' Error raading HEIGHTS . DAT I ' MARS116 

Purpose: 

To indicate a latitude read from HEIGHTS.DAT file was an unexpected value (file is 
missing, corrupted, or sized wrong based on parameters nlat and nlon). 

Remedy: 

Correct problem with HEIGHTS.DAT file. Rerun program. 


Code Producing the Message: 


260 Writ# (luO) 270 

270 formatt' Unabla to open Trajactory Data filat') 


MARS 14 4 
MARS145 


Purpose: 

To indicate an error occurred in attempt to open trajectory input data file 
TRAJDATA (file may not exist or may have wrong path set to its location). 

Remedy: 

Correct problem with TRAJDATA file. Rerun program. 
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Code Producing the Message: 

If (alsO . fft . 0. )Writa(iuO, *) • La out si da rang*. MARS177 

6 ' No dust storm assumed. ' HARS178 

Purpose: 

To indicate user attempted to simulate a dust storm for a value of areocentric 
longitude of the Sun (Ls) outside range 1 80 to 320. Message reminds user that dust storms 
cannot be simulated for Ls values outside this range and no dust storm is being simulated 
during the run. 

Remedy: 

None necessarily required. If a dust storm case is desired, select a different Ls (180 < 
Ls < 320). Rerun program. 


Code Producing the Message: 

WriteduO, *) ' Intensity must be between 0 and 3' MARS187 

Purpose: 

To indiciate a dust storm intensity < 0.0 or > 3.0 was entered. 

Remedy: 

Input a dust storm intensity between 0.0 and 3.0. Rerun program. 


Code Producing the Message: 

Write (iuO, *) ' F10.7 must be between 50 and 450' 


MARS213 


Purpose: 

To indicate an FI 0.7 solar flux <50 or >450 was entered. 
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Remedy: 


Input an FI 0.7 value between 50 and 450. Rerun program. 


Code Producing the Message: 

Write(iu0, *) ' Std. deviations must be between -3 and +3' MARS217 


Purpose: 

To indicate a standard deviation in the Stewart thermosphere model < -3 or > +3 
was entered (nominal value is 0.0). 

Remedy: 

Input a standard deviation value between -3 and +3. Rerun program. 


Code Producing the Message: 

430 Write (iuO, 440) 

440 FORMAT { ' Computing data. 1 ) 


MARS312 

MARS313 


Purpose: 

To indicate all data are input successfully and computations underway. 

Remedy: 

Complete the run. Time required depends on size of trajectory input file or number 
of automatic profile positions to be computed. 
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Code Producing the Message: 

9998 Stop ■ Error termination on opening or reading input file' MAKS 3 2 3 

Purpose: 

To indicate a termination message if file open error occurred on trajectory file or 
other files in FILES array. 

Remedy: 

Correct file problem, according to preceding error message. Rerun program. 


Code Producing the Message: 

9999 STOP ' Normal Termination' 

Purpose: 

To indicate successful completion of program run. 
Remedy: 

None required. 
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5.2 Batch Form Main Program MARSGRMB 


The batch form Marsgrmb main program does not produce diagnostic or progress 
messages other than the normal termination message for output to the STOP messages 
location. 


Code Producing the Message: 

999 STOP ' Normal Termination' 


MARB117 


Purpose: 

To indicate successful completion of program run. 
Remedy: 

None required. 
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5.3 Subroutine SETUP 


Several diagnostic or progress messages are produced by the SETUP subroutine. 
SETUP is called by the batch form only (not the interactive form) of the program. 


Code Producing the Message: 

50 Write (iuO, 60) files ( j ) , ierr ( j ) 

60 Format ( lx, al2, • File open error! Error =’,i5) 

Purpose: 

To indicate an error occurred in attempt to open one of the 17 files in FILES array 
(see separate Table 4-7 for file names). Message gives file name and error code that 
caused the (first encountered) problem. 

Remedy: 

Depending on meaning of error code on compiler system, take appropriate 
corrective action to ensure proper file open operation. Rerun program. 


SETU141 

SETU142 


Code Producing the Message: 


65 Stop ' Error or EOF on COSPAR.DAT file! ' 


SETU146g 


Purpose: 

To indicate an End-of-File or read error occurred on COSPAR.DAT file (file may 
not exist or path to this file may not be set appropriately). 

Remedy: 

Correct error condition (i.e., copy or restore file to right location). Rerun program. 
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Code Producing the Message: 

80 If (lat .ne. -100+10*i)Stop • Error reading HEIGHTS. DATl ' SETU150 

Purpose: 

To indicate a latitude read from HEIGHTS.DAT file was an unexpected value (file 
is missing, corrupted, or sized wrong based on parameters nlat and nlon). 

Remedy: 

Correct problem with HEIGHTS.DAT file. Rerun program. 


Code Producing the Message: 

Write ( iuO, 91) 

91 Format (' Input error in month, day or year.') 

Purpose: 

To indicate an illegal value of month, day, or year input via the NAMELIST 
INPUT file. 

Remedy: 

Correct erroneous value(s). Rerun program. 


SETU157 

SETU158 


Code Producing the Message: 


Write (iuO, 92) 

92 Format ( • Input error in hour, minute or seconds.' ) 


SETU170 
SETU17 1 


Purpose: 

To indicate an illegal value of hour, minute, or second input via the NAMELIST 
INPUT file. 


Remedy: 

Correct erroneous value(s). Rerun program. 



Code Producing the Message: 


260 

270 


Write (iuO, 270) 

format (' Unable to open Trajectory Data file! 1 ) 


SETU191 

SETU192 


Purpose: 


To indicate an error occurred in attempt to open trajectory input data file 
TRAJDATA (file may not exist or may have wrong path set to its location). 

Remedy: 


Correct problem with TRAJDATA file. Rerun program. 


Code Producing the Message: 

(alsO . gt . 0 .) Write ( iuO , # ) ' La outside range. 1 , 
& ' No dust storm assumed. 1 

Purpose: 


SETU216 

SETU217 


To indicate user attempted to simulate a dust storm for a value of areocentric 
longitude of the Sun (Ls) outside range 180 to 320. Message reminds user that dust storms 
cannot be simulated for Ls values outside this range and that no dust storm is being 
simulated during the run. ° 

Remedy: 

None necessarily required. If a dust storm case is desired, select a different Ls (180 
< Ls < 320). Rerun program. 


Code Producing the Message: 

Write (iuO, *) 1 Intensity must be between 0 3 1 


SETU222 


Purpose: 

To indicate a dust storm intensity < 0.0 or > 3.0 was entered. 
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Remedy: 


Input a dust storm intensity between 0.0 and 3.0 . Rerun program. 


Code Producing the Message: 

WritatiuO, *) ' riO. 7 must ba batwaan 50 and 450' SETU233 

Purpose: 

To indicate an FI 0.7 solar flux <50 or >450 was entered. 

Remedy: 

Input an FI 0.7 value between 50 and 450. Rerun program. 


Code Producing the Message: 

Writ* ( iuO , * ) ' Std. deviation* must b# batwaan -3 and +3' SETU237 


Purpose: 

To indicate a standard deviation in the Stewart thermosphere model <-3 or >+3 
was entered (nominal value is 0.0). 

Remedy: 

Input a standard deviation value between -3 and +3. Rerun program. 


Code Producing the Message: 

Writa(iu0, 291) 

291 Formate Error in parturbation modal numbar.') 


SZTU241 

SETU242 


Purpose: 

To indicate a value of the perturbation model code less than 1 or greater than 3 was 
input via NAMELIST INPUT file. 
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Remedy: 


Input a legal model number less than 1 or greater than 3. Rerun program. 


Code Producing the Message: 

Write (iuO, 292) 

292 Format ( * Error in starting random number. 1 ) 

Purpose: 


SETU247 

SETU248 


To indicate a random number seed value of less than 1 or greater than 29999 was 
input via NAMELIST INPUT file. 

Remedy: 


Input a random number seed value greater than 0 and less than 30000. Rerun 
program. 


Code Producing the Message: 

Write (iuO, 381) SETU295 

381 Format (■ x-code or y-code input error.') SETU296 

Purpose: 

To indicate a value of x-code outside the range 1 to 8, or y-code outside range 0 to 
8 was input via NAMELIST INPUT file. 

Remedy: 


Select x-code and y-code values for plotable output, according to the following 

codes: 


& 

& 

& 

& 

& 

& 

& 

& 

& 

& 

& 


Cod© Parameter'/ 

t/ 

1 Height (above reference ellipsoid, km)'/ 

2 Height (above local terrain, km) ' / 

3 Latitude (deg.)'/ 

4 West Longitude (deg.)'/ 

5 Time from start (Earth seconds)'/ 

6 Time from start (Martian Sols)'/ 

7 Areocentric Longitude of Sun, Ls (deg.)'/ 

8 Hour Angle for Local Time (Mars hours * 15)'// 

Use y-code = 0 for plotable output vs x-code variable only') 


SETU283 

SETU284 

SETU285 

SETU286 

SETU287 

SETU288 

SETU289 

SETU290 

SETU291 

SETU292 

SETU293 
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Rerun program. 


Code Producing the Message: 

Write (iuO, 382) SETU307 

382 Format (' Error in first latitude or longitude.') SETTJ308 

Purpose: 

To indicate a starting latitude absolute value greater than 90, or a starting longitude 
less than 0 or greater than 360, was input via NAMELIST input file. 

Remedy: 

Input a legal value of starting latitude and longitude. Rerun program. 


Code Producing the Message: 

9998 Stop • Error termination Check the LIST file for messages.' SETU32 

Purpose: 

To indicate termination, i.e., a fatal error condition occurred. Normally error 
conditions are described by message preceding this on screen output. For some options 
and error conditions, additional error information may appear in LIST output file. 

Remedy: 

Correct the error(s) noted. Rerun program. 


5-12 



o n 


5.4 Subroutines in the File MARSSUBS 


The subroutines in the MARSSUBS file produce diagnostic messages when certain 
error conditions occur. As an option, various diagnostic outputs are triggered by the user 
from within the Stewart thermosphere model, by setting FLAG (an integer) to a non-zero. 


Code Producing the Message: 

IF ( I FAULT .EQ. l)STOP ' PPND ERROR’ DSTP 61 

Purpose: 

To indicate a fatal error occurred in random number routine PPND. This is 
abnormal and may indicate a system-level problem. 

Remedy: 

Rerun the program with a different random number seed (NR1). If problem 
persists, run system level diagnostic tests to see if RANDOM and PPND perform properly. 


Code Producing the Message: 


9998 Stop ' Error termination reading data file I 1 


DSTP353 


Purpose: 

To indicate an error occurred while reading trajectory file TRAJDATA. 
Remedy: 

Correct data in TRAJDATA file. Rerun program. 


Code Producing the Message: 


if (FLAG . GT. 0 ) Write (iuO, 170) LSUN, DFAO, DZ 
170 FORMAT ( ' FROM PROC. DZDUST-- LSUN, DFAO, DZ « ’.3F8.2) 


DZDS 23 
DZDS 24 
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Purpose: 


This code is currently "commented out" (inactivated). If comment characters are 
removed from column 1 and if FLAG is set to a non-zero value and passed as an argument 
into subroutine DZDUST, then this optional diagnostic output is produced. 


Remedy: 

After tests are completed and diagnostic output is not desired, set FLAG to 0 (in 
Stewart2 subroutine, line STW2 20). 


Code Producing the Message: 


if (FLAO.gt.O)Writ#(iuO,*) ■ RREF, RAU, GZ - • ,RREF.RAU,OZ STW2 32 

if (FLAG.gt . 0)Writ#(iu0, * ) ' ZF,CHOT,ZZF « ■ , ZF,CHOT, ZZF STW2 44 

Writ#(iuO,*) ' RAU, FBARR, RREF, DR, DUST - ' STW2 49 

Writ#(iu0, 150)RAU,FBARR, RREF, DR, DUST STW2 50 

150 FORMAT (F6 . 2, F7 . 2, F8 .2, 2F6 . 2, F8 .2) STW2 51 

Writ#(iu0,*) 1 ' STW2 52 

Writ# ( iuO, *) 'FROM PROC. DRAG-- ZF,RF,TXNF,TF - ' STW2 53 

Writ#(iu0,160)ZF,RF,TINF,TF STW2 54 

160 FORMAT (F7 . 2 , F8 . 2 , 3F7 . 2 ) STW2 55 


Purpose: 

To produce this optional diagnostic output, set value of FLAG to a non-zero value. 


Remedy: 

After tests are completed and diagnostic output is not desired, set FLAG to 0 (in 
Stewart2 subroutine, line STW2 20). 


Code Producing the Message: 


if (FLAG .GT. 0) than 

Writa<iu0, 220 ) ZZF, TZ,MOLWTG, TOTALPRZ , TOTALNDZ , (NDZ(K) ,K-0, 3) , 
& TOTALMDZ, (MDZ(K) ,K=4,8) 

220 FORMAT ( 2F6 . 1, F5 . 1 , 6E10 . 3 / 17X, 6E10 . 3 ) 


THRM135 

THRM136 

THRM137 

THRM138 
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Purpose: 

To produce this optional diagnostic output, set the value of FLAG to a non-zero 
value and pass as an argument into subroutine THERMOS. 


Remedy: 

After tests are completed and diagnostic output is not desired, set FLAG to 0 (in 
Stewart2 subroutine, line STW2 20). 
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6. FUTURE PLANS 


The NASA Marshall Space Flight Center (MSFC) and NASA Jet Propulsion 
Laboratory (JPL) are working a joint program of future upgrades and applications 
development for Mars-GRAM. The following is a brief outline of the development 
features under consideration. 


6.1 Improved Mars Thermosphere Model 


The Stewart model thermosphere that serves as the basis for the upper altitudes in 
Mars-GRAM was designed to give global average values of density and other variables as 
a function of height and time. This means that latitude-longitude variations of 
thermospheric density, temperature, and pressure are unrealistic in Mars-GRAM. Since 
the wind components are based on extremely small horizontal pressure gradients (through 
the areostrophic or thermal wind relation), the winds estimated by Mars-GRAM in the 
thermosphere region are unrealistically small. 

To improve the Mars-GRAM thermospheric values, output data from the Bougher 
et al. (1988) Mars Thermospheric Global Circulation Model (MTGCM) will be used to 
derive appropriate modifier factors for latitude and time-of-day (longitude), for application 
to the Stewart model values. Comparisons with MTGCM results may also improve the 
Mars-GRAM (Stewart) model dependence on the effects of solar activity (10.7 cm solar 
flux). With comparison of wind output from MTGCM, the wind model in Mars-GRAM, 
including the viscous modification model, also has potential for improvement. 


6.2 Climatic Changes Since the Mariner-Viking (1970’s) Time Period 

Because Mars-GRAM was developed from parameterizations to atmospheric data 
observed by the Mariner and Viking missions, the model is representative of the Mars 
atmosphere during the 1970’s. Recent observations by Clancy et al. (1990) indicate, in 
response to atmospheric cooling as the dusty atmosphere has cleared in the last two 
decades, that current Mars temperature profiles are distinctly cooler than those observed in 
the Viking era. Comparisons of Mars-GRAM mid-latitude average temperature profiles 
with recent data from Clancy (provided by Rich Zurek, private communication), indicate 
about 20 K cooling at the 40-50 km height range, about 15 K at 20-30 km, but little change 
at 5-10 km. This temperature change could have significant effects on atmospheric density 
at high altitudes and, as Clancy et al. (1990) pointed out, is of considerable importance in 
planning for Mars missions that involve aerobraking. 



Consider Mars-GRAM profiles of temperature, density and pressure to be T 0 (z), 
Po(z), and po(z), respectively, with the perfect gas law and hydrostatic relations applicable, 
namely. 


Po = PoRT 0 


(23) 


and 


dpo/dz = -p 0 g = -p 0 (z)/H(z) . (24) 


Consider also that climatic changes up to the present time have induced new profiles of 
temperature, density and pressure, T|(z), p j (z) and pi(z). Application of equations (23) and 
(24) to both the original and the new atmospheric states yields the relation for density: 


Sp/Po = (P.-PoVPo = (To/Ti)cxpfj [(T, - To) / (T 0 H)]jdz' - 1 


(25) 


Consider the simple example in which temperature has uniformly decreased by 5% 
between the surface and 100 km, i.e., 8T / To = ( Tj - To ) / To = -0.05 at all altitudes. 
Since the scale height, H, is about 10 km, the integral term in equation (25) would have a 
value of about -0.5 and the density change at 100 km altitude from equation (25) would be 
5p / p 0 = -36%. In more general cases, the exact details of the density effect at altitude 
would, of course, depend on the profile of the temperature changes ( T, - T 0 ) and the 
profile of the scale height values H(z). 

Another Mars-GRAM enhancement under consideration by MSFC and JPL is to 
use new atmospheric temperature profiles from the Clancy data (and/or or from prior 
missions, after the first new Mars mission) to derive a climate-shift option (i.e. to estimate 
current conditions with the climatic shift option on, or to estimate Viking era conditions 
with the new climate shift option off). 


Mars-GRAM as an Operational Tool for Aerobraking 

MSFC and JPL are also interested in developing operational tools, based on Mars- 
GRAM, to use during the aerobraking phase of future Mars missions. Two potential 
approaches are to: 
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(1) Use Mars-GRAM densities evaluated near the periapsis position and time to 
compare with densities derived from the "real-time" orbit/trajectory analysis. Keep a time- 
series history of a "figure of merit" to monitor the trends of necessary adjustments to 
convert Mars-GRAM derived densities to actual densities for the aerobraking period. The 
figure of merit may be simply a ratio of Mars-GRAM value to the real time value or, 
perhaps, the value of the standard deviation (parameter STDL in the Stewart thermosphere 
model) required to get Mars-GRAM output to agree with real-time derived density. Once 
the observed trend in the figure of merit is established, use, in conjunction with Mars- 
GRAM calculations, to make short-term forecasts of the atmospheric density to be 
encountered over the next several passes through periapsis. 

(2) Use other "real-time" measurements of temperature profiles (e.g., Clancy-type 
ground-based microwave observations or temperature profiler data from sensors onboard 
the Mars orbiter) to derive estimates of necessary corrections to be applied to Mars- 
GRAM output. In this mode, estimates of density deviations from Mars-GRAM are 
computed via equation (25), using estimates of observed temperature profile deviations 
from Mars-GRAM profiles. 
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APPENDIX A 


DESCRIPTION OF THE MARS-GRAM PROGRAM AND SUBROUTINES 


This Appendix gives descriptions of the interactive form (Mars-GRAM) main 
program, the batch form (MARSGRMB) main program, the subroutine SETUP (used only by 
the batch form), and each of the functions and subroutines in the file MARSSUBS (in 
alphabetical order by the line number code for the MARSSUBS routines). References to 
functions, subroutines, and variables used appear as either all capital letters or as upper and 
lower case. The FORTRAN compiler is assumed to be case-insensitive. 

Each description gives: (1) function or subroutine name, (2) four-character line number 
code, (3) brief description of purpose of subroutine, (4) subroutine(s) that call function or 
subroutine, (5) common blocks used, and (6) list and descriptions of input variables, output 
variables, and local variables (those not passed through common blocks). For a list of 
variables in common blocks, see Table 2-3. 

A technical discussion of methods used by each function or subroutine is also given. 
Reference is made to the actual line of program code, including line numbers, for easy cross- 
reference to the program code. 
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The Mars-GRAM Main Program (Interactive Form) 


Main Program: Mars-GRAM (version 3.34) 

Description: Interactive form of the Mars-GRAM program. 


Cod*: MARS 


Call*d By: N/A 


C ommon Blocks Used (Se* Table 2-3 for a list of variables in common blocks): 

COSPARNH , DATACOM, FILENAME, RANDCOM, TERHGT , THERM, WAVEDAT 


Input Variables (Read from interactive input) : 


Name 

Type 

Description 

ALSO 

REAL *4 

time (Ls value) for beginning of dust storm 

DC 

REAL *4 

COSPAR density array 

DELHGT 

REAL *4 

height increment (km) 

DELLAT 

REAL *4 

northward latitude increment (degrees) 

DELLON 

REAL *4 

westward longitude increment (degrees) 

DELTIME 

REAL *4 

time increment (seconds) 

DUSTLAT 

REAL *4 

latitude of local-scale dust storm 

DUSTLON 

REAL *4 

longitude of local-scale dust storm 

F107 

REAL* 4 

10.7 cm solar flux at Earth position (1AU) 

FHGT 

REAL *4 

first height (km) 

FLAT 

REAL *4 

first latitude (degrees) 

FLON 

REAL *4 

first longitude (degrees West) 

IHR 

INTEGER* 4 

hour of day (UTC or GMT) 

IMIN 

INTEGER* 4 

minute of hour 

INTENS 

REAL *4 

intensity level for dust storm (0.0 to 3.0) 

LAT 

INTEGER* 4 

latitude for reading HEIGHTS.DAT input 

LOG SCALE 

INTEGER* 4 

output units for pressure and density (0=MKS, 

1 = log-base-10 MKS , 2=%deviation from COSPAR) 

LSTFL 

CHAR* 12 

Name of LIST file 

MDAY 

INTEGER* 4 

day of month 

MODPERT 

INTEGER* 4 

perturbation model (l=random, 2=wave, 3=both) 

MONTH 

INTEGER* 4 

month of year 

MYEAR 

INTEGER* 4 

year (4 digit; can be 2-digit if 1970 to 2069) 

NPOS 

INTEGER* 4 

maximum number of automatically-generated positions 

NR1 

INTEGER* 4 

random number seed (1-30,000) 

NVARX 

INTEGER* 4 

x-coordinate plot variable (See MARS262-MARS274) 

NVARY 

INTEGER* 4 

y-coordinate plot variable (See MARS2 62 -MARS274 ) 

OUTFL 

CHAR* 12 

Name of OUTPUT file 

PC 

REAL *4 

COSPAR pressure array 

RADMAX 

REAL *4 

maximum radius of dust storm (km; global if 
radmax=0 or radmax>10 , 000 ) 

SEC 

REAL *4 

second of minute 

STDL 

REAL *4 

number of standard deviations for thermospheric 
variation (-3.0 to +3.0) 

TC 

REAL* 4 

COSPAR temperature array 

TH 

REAL* 4 

terrain height array, read from HEIGHTS.DAT 

ZC 

REAL *4 

COSPAR height array 


A-2 



Output Variables (Includes echo out of some input variables) : 


Type Description 


ALS 

REAL *4 

(Ls) areocentric longitude of sun 

ALSO 

REAL *4 

time (Ls value) for beginning of dust storm 

DATE 

REAL* 8 

Julian date 

DUSTLAT 

REAL *4 

latitude of local-scale dust storm 

DUSTLON 

REAL *4 

longitude of local-scale dust storm 

F107 

REAL *4 

10.7 cm solar flux at Earth position (1AU); 
program also writes out value at Mars position 

I HR 

INTEGER* 4 

hour of day (UTC or GMT) 

IMIN 

INTEGER* 4 

minute of hour 

INTENS 

REAL *4 

intensity level for dust storm 

MDAY 

INTEGER* 4 

day of month 

MODPERT 

INTEGER* 4 

perturbation model (l=random, 2=wave, 3=both) 

MONTH 

INTEGER *4 

month of year 

MYEAR 

INTEGER* 4 

year (4 digit; can be 2-digit if 1970 to 2069) 

NPOS 

INTEGER* 4 

maximum number of automatically-generated positions 

NR1 

INTEGER* 4 

random number seed 

RADMAX 

REAL *4 

maximum radius of dust storm 

SEC 

REAL *4 

second of minute 

STDL 

REAL *4 

number of standard deviations for thermospheric 
variation 

THGT 

REAL *4 

local elevation of surface (km) 


Local 

Name 

Variables (not 

Type 

passed through commons) : 

Description 

ALS 


REAL *4 

(Ls) areocentric longitude of the sun 

CHGT 


REAL *4 

current height 

CLAT 


REAL* 4 

current latitude 

CLON 


REAL *4 

current longitude 

CSEC 


REAL *4 

current time from start (seconds) 

DATE 


REAL* 8 

Julian date 

DATEO 


REAL* 8 

initial Julian data 

DELHGT 


REAL *4 

height increment 

DELLAT 


REAL* 4 

latitude increment 

DELLON 


REAL *4 

longitude increment 

DELTIME 

REAL *4 

time increment 

DENS 


REAL *4 

atmospheric density 

DENSHI 


REAL *4 

high (approx. + 1 sigma) density 

DENSLO 


REAL*4 

low (approx. - 1 sigma) density 

DENSP 


REAL* 4 

perturbation in density (% of mean) 

EOF 


INTEGER* 4 

End-of-file flag for trajectory file 

EWWIND 


REAL *4 

eastward vector wind component (m/s) 

FHGT 


REAL *4 

first height (km) 

FILES 


CHAR* 12 

array of file names (see list of values elsewhere) 

FLAT 


REAL *4 

first latitude 

FLON 


REAL *4 

first longitude 

GZERO 


REAL *4 

gravity acceleration at z=0 

I 


INTEGER* 4 

index variable 

I DAY 


INTEGER* 4 

day-of-year for beginning of each month 

I ERR 


INTEGER* 4 

file open error flag array (equivalent to IERR1-IERR17 ) 

IERRl 


INTEGER* 4 

file open error flag 

IERR10 


INTEGER* 4 

file open error flag 

IERRl 1 


INTEGER* 4 

file open error flag 

IERR12 


INTEGER* 4 

file open error flag 

IERRl 3 


INTEGER*4 

file open error flag 

IERR14 


INTEGER* 4 

file open error flag 
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IERR15 

INTEGER* 4 

file open error flag 

IERR16 

INTEGER* 4 

file open error flag 

IERR17 

INTEGER* 4 

file open error flag 

IERR2 

INTEGER* 4 

file open error flag 

IERR3 

INTEGER* 4 

file open error flag 

IERR4 

INTEGER* 4 

file open error flag 

IERR5 

INTEGER* 4 

file open error flag 

IERR6 

INTEGER* 4 

file open error flag 

IERR7 

INTEGER* 4 

file open error flag 

IERR8 

INTEGER* 4 

file open error flag 

IERR9 

INTEGER* 4 

file open error flag 

I HR 

INTEGER* 4 

hour of day- 

IMIN 

INTEGER* 4 

minute of hour 

J 

INTEGER* 4 

index value 

L 

INTEGER* 4 

index value 

LAT 

INTEGER* 4 

latitude for reading terrain heights 

MARSAU 

REAL*4 

orbital position of Mars (in AU) 

MAXNUM 

INTEGER* 4 

maximum number of positions to compute 

MDAY 

INTEGER* 4 

day of month 

MONTH 

INTEGER* 4 

month of year 

MYEAR 

INTEGER* 4 

year 

NDAY 

INTEGER* 4 

day of year 

NLAT 

INTEGER* 4 

number of latitudes in terrain height array 

NLON 

INTEGER* 4 

number of longitudes in terrain height array 

NRl 

INTEGER* 4 

random number seed 

NSWIND 

REAL*4 

northward vector wind component 

PRES 

REAL*4 

atmospheric pressure 

RHO 

REAL* 4 

random number for sequence of perturbations 

SEC 

REAL*4 

seconds of minute 

SUNLAT 

REAL*4 

latitude of sub-solar position on surface 

SUNLON 

REAL*4 

longitude of sub-solar position on surface 

TEMP 

REAL*4 

atmospheric temperature (K) 

THGT 

REAL*4 

height of terrain surface above reference ellipsoid 

XYEAR 

REAL* 8 

intermediate year for Julian day calculation 


Methodology : 


The program prompts the user for values of the necessary input . 
Depending on the options selected, some input parameters assume default 
values and it is not necessary to input a value. Several input variables 
are also echoed back as output variable on the LIST file. All subroutines 
required by Mars-GRAM interactive form are contained in the marssubs file. 

Mars-GRAM opens all input and output files, interactively reads the 
necessary input values and steps through the output positions, generating 
and writing the output. If the trajectory option is selected (NPOS=0) , 
positions are read from the TRAJDATA file, until and end-of-file is 
encountered. If an automatically-generated profile of positions is used 
(NPOS>0) , output continues until the maximum allowed number of output 
positions (MAXNUM=NP0S-1) has been exceeded. 

If the time is within the period for dust storm formation (180 < Ls 
<320, Ls = areocentric longitude of sun), then a local-scale or a global- 
scale dust storm can be simulated. Position and or size of the dust storm 
are controlled by input values of DUSTLAT, DUSTLON, RADMAX and INTENS. 
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The Mars-GRAM Main Program (Batch Form) 


Main Program: Marsgrmb (version 3.34) Code: MARB 

Description: Batch form of the Mars-GRAM program. 

Called By: N/A 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

None 


Input Variables : 

There are no input variables directly to Marsgrmb. Input data are read, in 
NAMELIST form, by the subroutine SETUP. An example of the NAMELIST input 
is as follows: 


$ INPUT 
LSTFL 
OUTFL 
MONTH 
MDAY 
MYEAR 
NPOS 

I HR 
IMIN 
SEC 
ALSO 

INTENS 

RADMAX 

DUSTLAT 

DUSTLON 

F107 

STDL 

MODPERT 

NR1 

NVARX 

NVARY 

LOG SCALE 

FLAT 

FLON 

FHGT 

DELHGT 

DELLAT 

DELLON 

DELTIME 

$END 


= 'LIST', 

= 'OUTPUT' 
= 7, 

= 20 , 

= 76, 

= 21 , 

= 12 , 

= 30, 

= 0 . 0 , 

= 0 . 0 , 

= 0 . 0 , 

= 0 . 0 , 

= 0 . 0 , 

= 0 . 0 , 

= 185.0, 

= 0 . 0 , 

= 3, 

= 1001 , 

= 2 , 

= 0, 

= 0, 

= 22 . 0 , 

- 48.0, 

= -0.5, 

= 10 . 0 , 

= 0 . 0 , 

= 0 . 0 , 

= 0 . 0 , 


List file name (CON for console listing) 
Output file name 
month of year 
day of month 

year (4-digit; 1970-2069 can be 2-digit) 
max # positions to evaluate (0 = read data 

from TRAJDATA file) 

GMT hour of day 
minute of hour 

second of minute (for initial position) 
starting Ls value (degrees) for dust storm 

(0 = none) 

dust storm intensity (0.0 - 3.0) 
max. radius (km) of dust storm (0 or 

>10000 = global) 

latitude (deg) for center of dust storm 
West longitude (deg) of center of dust storm 
10.7 cm solar flux (10**-22 W/cm**2 at 1 AU) 
std. dev. for thermosphere variation (-3.0 

to +3.0) 

perturbation model; l=random, 2=wave, 3=both 
starting random number (0 < NRl < 30000) 
x-code for plotable output (l=hgt above ref. 

ellipse) 

y-code for 2-D plotable output (0 for 1-D 

plots ) 

1 for log-base-10 scale plots, 0 for linear 
scale, 2 for % deviation from COSPAR 
initial latitude (N positive) , degrees 
initial longitude (West positive) , degrees 
initial height (km) , above ref. ellipse 
height increment (km) between steps 
latitude increment (deg) between steps 
West longitude increment (deg) between steps 
time increment (sec) between steps 


MARB 18 
MARB 19 
MARB 20 
MARB 21 
MARB 22 
MARB 23 
MARB 24 
MARB 25 
MARB 26 
MARB 27 
MARB 28 
MARB 29 
MARB 30 
MARB 31 
MARB 32 
MARB 33 
MARB 34 
MARB 35 
MARB 36 
MARB 37 
MARB 38 
MARB 39 
MARB 40 
MARB 41 
MARB 42 
MARB 43 
MARB 44 
MARB 45 
MARB 46 
MARB 47 
MARB 48 
MARB 49 
MARB 50 
MARB 51 
MARB 52 
MARB 53 


Output Variables: 

There are no output variables that are produced directly by Marsgrmb. 
Standard LIST and OUTPUT files (as well as the output files for graphics 
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input, units 21-34) may be produced by subroutine DATASTEP. The Marsgrmb 
batch form is designed so that the output variables from DATASTEP may 
easily be output by the user (or passed to other user subroutines), in 
whatever format is desired. Output variables from DATASTEP are: 


TEMP = temperature (K) MARB 98 
PRES = pressure (N/m**2) MARB 99 
DENSLO = nominal low density (kg/m**3), approx. -1 sigma MARB100 
DENS = mean density (kg/m**3) MARB101 
DENSHI = nominal high density (kg/m**3), approx. +1 sigma MARB102 
DENSP = density perturbation about mean (% of mean) MARB103 
EWWIND = eastward wind component (m/s) MARB104 
NSWIND = northward wind component (m/s) MARB105 


If none of the LIST file or OUTPUT file or graphics file output is desired, 
this is accomplished by setting file unit parameter iup = 0, in line BLKD 
16. 


Local Variables (not passad through commons) : 


Name 

Type 

Description 

CHGT 

REAL *4 

current height 

CLAT 

REALM 

current latitude 

CLON 

REAL* 4 

current longitude 

CSEC 

REALM 

current time from start (seconds) 

DATEO 

REALM 

initial Julian date 

DELHGT 

REAL *4 

height increment 

DELLAT 

REALM 

latitude increment 

DELLON 

REAL* 4 

longitude increment 

DELTIME 

REAL* 4 

time increment 

DENS 

REALM 

atmospheric density 

DENSHI 

REALM 

high (approx. + 1 sigma) density 

DENSLO 

REAL *4 

low (approx - 1 sigma) density 

DENSP 

REAL *4 

perturbation in density (% of mean) 

EOF 

INTEGER* 4 

End-of-file flag for trajectory file 

EWWIND 

REALM 

eastward vector wind component 

I 

INTEGER* 4 

index variable 

MAXNUM 

INTEGER* 4 

maximum number of positions to compute 

NSWIND 

REALM 

northward vector wind component 

PRES 

REALM 

atmospheric pressure 

RHO 

REALM 

random number for sequence of perturbations 

TEMP 

REALM 

atmospheric temperature (K) 


Methodology: 

Stripped of all of the embedded comment statements, the batch form 
Marsgrmb consists of only the following lines of code. All variables are 
initialized through the SETUP subroutine and parameters are evaluated at 
successive locations (either trajectory or automatic profile mode) with the 
DATASTEP subroutine. Any desired output (or transfer to other subroutines) 
can be done within the loop that calls DATASTEP. 

DOUBLE PRECISION DATEO 
Real NSWIND 
Integer EOF 

Call Setup ( CHGT , CLAT , CLON , CSEC , DATEO , RHO , DELHGT , DELLAT , DELLON , 

& DEL TIME , MAXNUM) 

DO 900 I = 0, MAXNUM 

Call Datastep ( I , CHGT, CLAT , CLON, CSEC , DATEO , RHO, EOF , DELHGT , 

& DELLAT , DELLON , DELTIME , TEMP , PRES , DENSLO , DENS , DENSHI , DENSP , 

& EWWIND, NSWIND) 
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If (EOF .eq. 1 ) Goto 999 
900 Continue 

999 STOP ' Normal Termination' 
END 
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3 . 


Description of the SETUP Subroutine (for the Batch Form) 


Subroutine : SETUP Code: SETU 

Description: Opens files and reads input data for the Marsgrmb batch form 
program . 

Called By: Marsgrmb (main) 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

COSPARNH , DATACOM, FILENAME , RANDCOM, TERHGT, THERM 

Input Variables : 

Input data are read, in NAMELIST form, by the subroutine SETUP. An 
example of the NAMELIST input is as follows: 


$ INPUT 




MARB 

18 

LSTFL 

= 

’LIST' , 

! List file name (CON for console listing) 

MARB 

19 

OUTFL 

zr 

1 OUTPUT * 

! Output file name 

MARB 

20 

MONTH 

= 

7, 

! month of year 

MARB 

21 

MDAY 

= 

20, 

! day of month 

MARB 

22 

MYEAR 

- 

76, 

! year (4-digit; 1970-2069 can be 2-digit) 

MARB 

23 

NPOS 

- 

11, 

! max # positions to evaluate (0 = read data 

MARB 

24 




! from TRAJDATA file) 

MARB 

25 

IHR 

- 

12, 

! GMT hour of day 

MARB 

26 

IMIN 

— 

30, 

! minute of hour 

MARB 

27 

SEC 

= 

0.0, 

! second of minute (for initial position) 

MARB 

28 

ALSO 

— 

0.0, 

! starting Ls value (degrees) for dust storm 

MARB 

29 



i (0 = none) 

MARB 

30 

INTENS 

= 

0.0, 

! dust storm intensity (0.0 - 3.0) 

MARB 

31 

RADMAX 

= 

0.0, 

! max. radius (km) of dust storm (0 or 

MARB 

32 


! >10000 = global) 

MARB 

33 

DUSTLAT 

— 

0.0, 

! latitude (deg) for center of dust storm 

MARB 

34 

DUSTLON 

= 

0.0, 

! West longitude (deg) of center of dust storm 

MARB 

35 

F107 

- 

185.0, 

! 10.7 cm solar flux (10**-22 W/cm**2 at 1 AU) 

MARB 

36 

STDL 

- 

0.0, 

! std. dev. for thermosphere variation (-3.0 

MARB 

37 


i to +3 . 0) 

MARB 

38 

MODPERT 

— 

3, 

! perturbation model; l=random, 2=wave, 3=both 

MARB 

39 

NR1 

- 

1001, 

! starting random number (0 < NRl < 30000) 

MARB 

40 

NVARX 

_ 

1, 

! x-code for plotable output (l=hgt above ref. 

MARB 

41 


1 ellipse) 

MARB 

42 

NVARY 

- 

o, 

t y-code for 2-D plotable output (0 for 1-D 

MARB 

43 


j plots) 

MARB 

44 

LOGSCALE 

= 

o, 

! 1 for log-base-10 scale plots, 0 for linear 

MARB 

45 



! scale, 2 for % deviation from COSPAR 

MARB 

46 

FLAT 

= 

22.0, 

! initial latitude (N positive) , degrees 

MARB 

47 

FLON 

- 

48.0, 

j initial longitude (West positive) , degrees 

MARB 

48 

FHGT 


-0.5, 

! initial height (km), above ref. ellipse 

MARB 

49 

DELHGT 

= 

10.0, 

! height increment (km) between steps 

MARB 

50 

DELLAT 

= 

0.0, 

! latitude increment (deg) between steps 

MARB 

51 

DELLON 

— 

0.0, 

! West longitude increment (deg) between steps 

MARB 

52 

DELTIME 

= 

0.0, 

! time increment (sec) between steps 

MARB 

53 


$END 
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Output Variables (not passed through commons): 


Name 

Type 

Description 

CHGT 

REAL *4 

current height 

CLAT 

REAL *4 

current latitude 

CLON 

REAL *4 

current longitude 

CSEC 

REAL* 4 

current time (seconds) 

DATEO 

REAL* 4 

initial Julian date 

DHGT 

REAL*4 

height increment 

DLAT 

REAL* 4 

latitude increment 

DLON 

REAL*4 

longitude increment 

DTI ME 

REAL*4 

time increment (seconds) 

MAXNUM 

INTEGER* 4 

maximum number of positions to compute 

RHO 

REAL* 4 

random number for sequence of perturbations 


Local Variables (not passad through commons) : 


Name 

Type 

Description 

ALS 

REALM 

areocentric longitude of the sun (Ls) 

DATE 

REAL* 8 

Julian date 

DELHGT 

REAL *4 

height increment 

DELLAT 

REAL* 4 

latitude increment 

DELLON 

REAL* 4 

longitude increment 

DELTIME 

REAL *4 

time increment 

FILES 

CHAR* 12 

file name array (see Table with list of files) 

GZERO 

REAL*4 

acceleration of gravity at z = 0 

I 

INTEGER* 4 

index variable 

IDAY 

INTEGER* 4 

day-of-year for beginning of each month 

I ERR 

INTEGER* 4 

file open error flag array (equivalent to IERR1-IERR17) 

IERRl 

INTEGER* 4 

file open error flag 

IERR10 

INTEGER* 4 

file open error flag 

IERRl 1 

INTEGER* 4 

file open error flag 

IERR12 

INTEGER* 4 

file open error flag 

IERR13 

INTEGER* 4 

file open error flag 

IERR14 

INTEGER* 4 

file open error flag 

IERRl 5 

INTEGER* 4 

file open error flag 

IERRl 6 

INTEGER* 4 

file open error flag 

IERR17 

INTEGER* 4 

file open error flag 

IERR2 

INTEGER* 4 

file open error flag 

I ERR 3 

INTEGER* 4 

file open error flag 

IERR4 

INTEGER* 4 

file open error flag 

IERR5 

INTEGER* 4 

file open error flag 

I ERR 6 

INTEGER* 4 

file open error flag 

IERR7 

INTEGER* 4 

file open error flag 

IERR8 

INTEGER* 4 

file open error flag 

IERR9 

INTEGER* 4 

file open error flag 

J 

INTEGER* 4 

index variable 

L 

INTEGER* 4 

index variable 

LAT 

INTEGER* 4 

latitude for reading terrain heights 

MARSAU 

REAL *4 

orbital position of Mars (in AU) 

NDAY 

INTEGER* 4 

day of month 

NLAT 

INTEGER* 4 

number of latitudes in terrain height array 

NLON 

INTEGER* 4 

number of longitudes in terrain height array 

SUNLAT 

REAL *4 

latitude of sub-solar position on surface 

SUNLON 

REAL *4 

longitude of sub-solar position on surface 

THGT 

REAL *4 

height of terrain surface above reference ellipsoid 

XYEAR 

REAL *8 

intermediate year for Julian day calculation 
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Methodology: 

Reads NAMELIST input values (see above) , computes Julian day, open output 
files and passes data to the Marsgrmb main program through the common blocks . 
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4. Description of the MARSSUBS.FOR Functions and Subroutines 


Function: 


Alb 


Code s ALBL 


Description ; Surface albedo as a function of areocentric longitude of the 
sun (Ls) and latitude. 

Called By: Tsurface 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

None 


Input Variables (not passed through commons) : 


Name 


Type 


Description 


als . . . REAL *4 

alat. . . REAL* 4 


areocentric longitude of sun (Ls), degrees 
latitude, degrees 


Output Variables (not passed through commons) : 

Name Type Description 

Alb . . . REAL*4 surface albedo 


Local Variables (not passed through commons): 


Name 


Type Description 


al . . 

also . 
caprad 
cldrad 
pil80 
radc . 


REAL *4 intermediate variable for albedo calculation 

REAL*4 Ls phase for polar hood cloud 

REAL* 4 polar cap radius, degrees 

REAL*4 polar hood cloud radius, degrees 

REAL*4 factor for degrees to radians 

REAL*4 radius amplitude for polar hood cloud 


Methodology : 

Computes albedo of non-polar surface as 


al = 0.32 - 0.12* cos (pil80*alat) 


ALBL 9 


Computes radius for polar hood cloud from 


cldrad = radcMl. + cos (pil80* (als - alsc))) 


ALBL 22 


where alsc = 280., radc = 20. for North polar hood cloud and alsc = 60. and radc = 
27. for South polar hood cloud. If the latitude is poleward of the polar hood cloud 
radius, the surface albedo is increased by 15 % 


If (abs (alat) .ge. 90. - cldrad)al = 1.15*al 


ALBL 24 


A-ll 


The radius of the polar cap is computed from the polecap function 


caprad = polecap (alat , als ) 


ALBL 26 


If the latitude is poleward of the polar cap radius, the surface albedo is set to 

0.6 

If (abs(alat) .ge. 90. - caprad) al = 0.6 ALBL 28 
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Function: ampint Cods: 

®** cr iption: Interpolate wave amplitudes, linearly on input variable 
Called By: Wavepert 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

None 

Input Variables (not passed through commons) : 

Name Type Description 

ampl . . . REALM amplitude at 1st value of input variable (xl) 

amp2 . . . REALM amplitude at 2nd value of input variable (x2) 

dx - • • • REAL *4 relative difference in input [ (x-xl) / (x2-xl) ] 

Output Variables (not passed through commons) : 

Name Type Description 

ampint. . REALM interpolated value of wave amplitude 

Local Variables (not passed through commons) : 

Name Type Description 

None 

Methodology : 

Interpolates linearly to position x, between ampl (at xl) and amp2 (at x2 ) by 
the relation 

ampint = ampl + (amp2 - ampl)*dx AMPN 2 

where dx = (x - xl)/ (x2-xl) . The input variable x may be height or latitude. 
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Subroutine: ATMOS2 


Code: ATM2 


Description: Driver routine to call either the STEWART 2 thermosphere model 
or the stratos, lower altitude model, depending on the current height 

Called By: Datastep 

C ommo n Blocks Ussd (See Table 2-3 for a list of variables in c ommo n blocks) : 

THERM 


Input Variables (not passed through commons): 


Name 

Type 

Description 

ALS. . . 

REALM 

areocentric longitude of the sun (Ls) , degrees 

ALSO . . 

REAL* 4 

initial value of areocentric longitude of the sun 

CHGT . . 

REAL* 4 

current height (km) 

CLAT . . 

REALM 

current latitude (degrees) 

CLON . . 

REAL *4 

current West longitude (degrees) 

DUSTA . . 

REAL* 4 

dust storm magnitude (0-1) for daily amplitude 

DUSTM . . 

REAL* 4 

dust storm magnitude (0-1) for average T and p 

INTENS . 

REAL* 4 

dust storm intensity (0.0 - 3.0) 

IUO . . . 

INTEGER* 4 

unit number for screen output messages 

MARSAU . 

REAL* 4 

Mars orbital radius (AU) 

SUNLAT . 

REAL *4 

latitude of sub-solar position (degrees) 

SUNLON . 

REAL *4 

West longitude of sub-solar position (degrees) 

20 ... . 

REAL*4 

local terrain height rel . to ref. ellipsoid (km) 


Output Variables (not passed through commons): 


Name 

Type 

Description 

BRUNTF . . 

REAL *4 

Brunt-Vaisala frequency 

DENST . . . 

REAL *4 

atmospheric density (kg/m^) at output position 

DENSURF . . 

REAL* 4 

atmospheric density (kg/m ) at surface 

H 

REAL *4 

scale height (km) at output position 

LWFCTR . . 

REALM 

lower deviation factor foj density 

PRES . . . 

REAL* 4 

atmospheric pressure (N/m ) 

RSC .... 

REAL* 4 

areocentric radius to output position (km) 

TAVG . . . 

REAL *4 

daily average surface temperature (K) 

TEMP . . - 

REALM 

temperature at output position (K) 

TFACTOR . . 

REAL* 4 

density perturbation factor at base of thermosphere 

TMAX . . . 

REAL* 4 

daily maximum surface temperature (K) 

TMIN . . . 

REAL *4 

daily minimum surface temperature (K) 

UPFCTR . . 

REALM 

upper deviation factor for density 

ZF .... 

REALM 

height of base of thermosphere (km) 


Local Variables (not passed through commons): 


Name Type 


Description 


AMF . 
AMHI . 
AML0. 
AMW . 
BVF . 
BVFO 
DENSHI 
DENSLO 


REAL* 4 
REAL* 4 
REAL* 4 
REAL *4 
REAL *4 
REAL*4 
REAL* 4 
REAL *4 


molecular weight at base of thermosphere 

molecular weight for high density 

molecular weight for low density 

average molecular weight, surface to turbopause 

Brunt-Vaisala frequency 

Brunt-Vaisala frequency at surface 

high (approx. +1 sigma) density 

low (approx. -1 sigma) density 


A- 14 



DLON. . 


REAL *4 

DR. . . 


REAL *4 

DUST . . 


REAL *4 

ES. . . 


REAL* 4 

GAM . . 


REAL *4 

GAM0 . . 


REAL *4 

GOR . . 


REAL *4 

GPH . . 


REAL *4 

GPHGT . 


REAL* 4 

GZ . 


REAL *4 

GZERO . 


REAL* 4 

HHI . . 


REAL *4 

HLO . . 


REAL* 4 

I . . . 


INTEGER* 4 

J . . . 


INTEGER* 4 

P . . . 


REAL *4 

PAVG . . 


REAL *4 

PF. . . 


REAL *4 

PFAC . . 


REAL *4 

PRESHI . 


REAL *4 

PRESLO . 


REAL *4 

PSURF . 


REAL *4 

RBAR . . 


REAL *4 

RREF . . 


REAL *4 

RSTAR . 


REAL* 4 

SHGT . . 


REAL *4 

SMA . . 


REAL *4 

T . . . 


REAL *4 

T0BAR . 


REAL *4 

TEMPHI . 


REAL *4 

TEMPLO . 


REAL *4 

TF. . . 


REAL *4 

TGRAD . 


REAL *4 

TIME. . 


REAL *4 

TO. . . 


REAL *4 

TSURF . 


REAL *4 

Z . . . 


REAL* 4 


longitude difference between sun and output position 

Correction to ZF for pressure variation 

Correction to ZF for dust storm 

long-term standard deviations about nominal ZF 

temperature lapse rates between significant heights 

mean lapse rates, for use in perturbation model 

gravity divided by gas constant 

areopotential heights at significant levels 

areopotential height difference for interpolation 

gravity at output position 

surface gravity 

scale height for high density 

scale height for low density 

index variable 

index variable 

array of pressures at significant levels 

daily average surface pressure 

intermediate value in pressure calculation 

Seasonal relative pressure variation 

pressure for high density 

pressure for low density 

local position surface pressure 

gas constant for Martian atmosphere up to turbopause 
erence radius to Cain 6 . 1 mb surface 
Universal gas constant 
height for call to thermosphere model 
semi-major axis of Mars orbit 
array of temperatures at significant levels 
average surface temperature without dust storm effect 
temperature for high density 
temperature for low density 

intermediate value in temperature calculation 
temperature gradient used in Brunt-Vaisala calculation 
local solar time 

factor used in correction to ZF for pressure variation 
average surface temperature wit dust storm effect 
areographic heights at significant levels 


Methodology ; 


Evaluates surface temperature (K) , by calling Tsurface (ATM2 56). Evaluates 
acceleration of gravity at altitude CHGT by calling RELLIPS (ATM2 59). Gets mean 
temperature lapse rates, for use in perturbation model, by calling Gamma (ATM2 62). 

Evaluates temperatures at significant levels (0.,5.,15 30 50 75 

areopotential km) by calling Temps { ATM2 69). Evaluates surface pressure / N/m 2 , by 
calling Psurface (ATM 2 71) . Evaluates pressure at significant levels by calling 
Pressure {ATM 2 74) . Gets surface density from perfect gas law 


densurf = psurf/ (Rbar*Tsurf ) 


ATM2 76 


Computes Brunt Vaisala frequency at surface 


bvfO = gzero* ( tgrad + gzero/Cp (Tsurf ) ) /Tsurf 
if (bvf 0 .It. 0.0) bvf 0 = 1.0E-5 
bvf 0 = Sqrt(bvfO) 


ATM2 79 
ATM 2 80 
ATM 2 81 


the current hei ^ ht is below -5 km, the surface conditions are returned and 
the subroutine is exited (ATM2 82 - ATM2 93). If the current height is below z(5) 
(areographic height- for 75 km areopotential height), ATMOS 2 uses the lapse rates 
from the significant levels to interpolate for temperature, pressure and density 
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o n 


TEMP = T(j) - gam( j+1) *gphgt 
If (abs (gam( j+1) ) .gt. 0.001)then 

PRES = p ( j ) *exp (goR*alog (TEMP/T ( j ) ) /gam( j+1) ) 
Else 

PRES = p ( j ) *exp (-goR*gphgt/T ( j ) ) 

Endif 

DENST = PRES/ (Rbar*TEMP) 


ATM2101 

ATM2102 

ATM2103 

ATM2104 

ATM2105 

ATM2106 

ATM2107 


and solves for the Brunt-Vaisala frequency and high and low density factors 


bvf = gzMtgrad + gz/Cp (TEMP) ) /TEMP 
if (bvf .le. 0.0)bvf = 1.0E-5 
Brunt f = Sqrt(bvf) 

UPFCTR = density and temperature part of un-damped mountain 
waves 

UPFCTR = (bvfO*Bruntf /gz) *Sqrt ( dens urf /DENST) 

LWFCTR =0.0 

If the current height is above z(5), then conditions are evaluated at height 
ZF, the base of the thermosphere (ATM2143 and ATM2146) . If the current height is 
above z(5) but below ZF, values are found by interpolation, using the Stratos 
subroutine (ATM2150 - ATM2163) . If the current height is above the base of the 
thermosphere, values are computed from the STEWART2 thermosphere model (ATM2164 - 
ATM2172 ) . 


ATM2112 

ATM2113 

ATM2114 

ATM2115 

ATM2116 

ATM2117 

ATM2118 
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Subroutine : Block Data 


Code: BLKD 


Description: Loads data values into common blocks, for use by various 
subroutines in the program 

Called By: N/A 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

DATACOM, WAVEDAT 

Input Variables (not passed through commons) : 

Name Type Description 

None 

Output Variables (not passed through commons) : 

Outputs the variables in the commons DATACOM and WAVEDAT 


Local Variables (not passed through commons): 

Name Type Description 

None 


Methodology : 

Uses data statements to fill the arrays and assign values to the 
variables in the commons DATACOM and WAVEDAT. These commons are used by 
several subroutines (see Table 2-2) . 
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Subroutine: cospar 


Code: COSP 


Description: Computes COSPAR northern hemisphere mean values of pressure, 
density and temperature as a function of height. 

Called By: Datastep 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

cosparnh 

Input Variables (not passed through commons) : 

Name Type Description 

z REAL *4 height (km) at which to evaluate COSPAR data 

Input arrays of COSPAR temperature (tc, K) , pressure (pc, mb) and 
density (dc, g/m**3) come from the common COSPARNH, from data in the 
CO SPAR . DAT file. 


Output Variables (not passed through commons): 


Name Type 


Description 


p REAL* 4 

rho. . . . REAL *4 

t REAL * 4 


COSPAR pressure (N/nt) at height z 
COSPAR density (kg/m 3 ) at height z 
COSPAR temperature (K) at height z 


Local Variables (not passed through commons): 


Name Type 


Description 


aexp 
dz . 
H. . 
iz . 
R. . 
R1 . 
R2 . 


REAL *4 
REAL *4 
REAL *4 
INTEGER* 4 
REAL *4 
REAL *4 
REAL *4 


exponent for constant lapse rate interpolation 
relative height displacement for interpolation 
scale height for isothermal interpolation 
height index for interpolation 
interpolated gas law constant 
gas constant at lower interpolation height 
gas constant at upper interpolation height 


Mathodology: 

Uses an array of COSPAR data values (loaded in from the COSPAR.DAT 
file) . The height interval is 1 km from -5 to 130 km and 10 km from 130 to 
360 km. Computed height index, iz, is for the COSPAR height value just 
below the desired height z (COSP 10 - COSP 14) . Return zero values if the 
height is outside the range -5 to 360 km. (COSP 16 - COSP 21) . 

Uses linear interpolation on temperature 


dz = (z - zc (iz) ) / (zc (iz+1) - zc(iz)) 

t = tc(iz) + (tc(iz+l) - tc(iz))*dz 


COSP 24 
COSP 25 


Uses a constant lapse rate, power-law interpolation (unless the layer is isothermal) 
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aexp = Alogfpc (iz+1) /pc(iz) ) /Alog ( tc ( iz+1 ) /tc(iz) ) 
p = 100. *pc<iz) * (t/tc(iz) ) **aexp 


Uses an exponential (constant scale height) interpolation if isothermal 
prevail 


H = (zc (iz+1) - zc (iz) ) /Alog (pc (iz) /pc (iz+1) ) 
p = 100 . *pc (iz) *exp (- (z-zc (iz) ) /H) 


Computes density by linearly interpolating on gas constant and applying 
gas law relation 


Rl = pc (iz) / (dc(iz) *tc(iz) ) 

R2 = pc (iz+1) / (dc (iz+1) *tc (iz+1) ) 

R = Rl + (R2-R1) *dz 

density from perfect gas law (and convert units to kg/m**3) 
rho = 10 . *p/ (R* t) 


COSP 28 
COSP 2 9 


conditions 


COSP 31 
COSP 32 


the perfect 


COSP 35 
COSP 36 
COSP 37 
COSP 38 
COSP 3 9 
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Functions Cp 


Code: CPOT 


Description: Specific heat at constant pressure, for a C02 atmosphere, as 
function of temperature. 

Called By: ATMOS 2 

Common Blocks Used (See Table 2-3 for a list of variables in c ommo n blocks): 

None 

Input Variables (not passed through commons) : 

Name Type Description 

T . . . . REAL *4 temperature (K) 

Output Variables (not passed through commons): 

Name Type Description 

Cp. . . . REAL *4 specific heat (J kg 1 K 1 ) 

Local Variables (not passed through commons) : 

Name Type Description 

None 

Methodology: 

Computes Cp{T) from formula 


Cp = 639.5 + 0 . 123687*T + 0 . 00200225*T*T 


CPOT 4 



Subroutine ; Du s t f ac t 


Coda : DSTF 


Dascription: Computes relative dust storm intensity factors dustM and dustA as a 
function of the time since start of the storm, (als - alsO) , measured in Ls angle 
(degrees), and as a function of the storm intensity, intens . dustM is for relative 
magnitude of effect on daily average temperature and pressure. dustA is for 
relative magnitude of effect on diurnal (and semi-diurnal) amplitudes for 
temperature and pressure. 

Called By: Datastep 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

DATACOM 

Input Variables (not passed through commons) : 


Name 

Type 

Description 


ALS . . 

. . REAL *4 

Areocentric longitude of sun (Ls) 


CLAT . . 

. . REAL* 4 

current latitude 


CLON . . 

. . REAL* 4 

current longitude 


HGT . . 

. . REAL* 4 

current height 


Output 

Variables (not passed through commons) : 


Name 

Type 

Description 


DUSTA . 

. . REAL *4 

relative magnitude for dust storm effects 
variation amplitudes for temperature and 

on daily 
1 pressure 

DUSTM . 

. . REAL* 4 

relative magnitude for dust storm effects 
average temperature and pressure 

on daily 

Local Variables (not 

passed through commons) : 


Name 

Type 

Description 



DEW . . 

DLS . . 

DNS . . 

HGTDUST 
RAD . . 

RADDUST 
SIZEFACT 


REAL *4 
REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 


Methodology : 


east-west component of distance from dust storm center 

Ls difference from starting Ls (LsO) for dust storm 

north-south distance component from dust storm center 

vertical size to top of dust storm 

radial distance from dust storm center 

full horizontal radius of dust storm size 

shape factor to compute DUSTA and DUSTM at location 


Computes Ls difference (dls, degrees) from Ls value at start of dust storm 
(LsO), and returns zero dust storm effect if dls < 0 or dls > 48. 


dls = als - alsO 

If (dls .le. 0.0 .or. dls .gt. 48.0)then 
dustM = 0.0 
dustA = 0.0 
Return 
End if 


DSTF 14 
DSTF 15 
DSTF 16 
DSTF 17 
DSTF 18 
DSTF 19 
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Computes an initial factor value dustM that increases linearly from 0 to 1 as dls 
increases from 0 to 6 degrees, decreases from 1 back to 0 as dls varies between 24 
and 48 degrees, or dustM = 1 otherwise. 

If (dls .le. 6.0) then 
dustM = dls/ 6. 

Else If (dls .ge. 24.) then 
dustM = 2. - dls/24. 

Else 

dustM = 1.0 
Endif 

Computes an initial factor value dustA that increases linearly from 0 to 1 as dls 
increases from 0 to 9, decreases from 1 back to 0 as dls varies between 18 and 48 
degrees, or dustA = 1 otherwise. 

If (dls .le. 9.0) then 
dustA = dls/ 9 . 0 
Else If (dls .ge. 18 .) then 
dustA = (48. - dls)/30. 

Else 

dustA = 1.0 
Endif 

Computes the size factor, sizefact, based on the position from the center of the 
dust storm and the maximum horizontal and vertical size of the storm. 

sizefact = 1.0 
If (radmax .ne. 0.0) Then 
sizefact = 0.0 

dns = DTR*Rref * (CLAT - dustlat) 
dew = DTR*Rref *cos (DTR*CLAT) * (CLON - dustlon) 
rad = Sqrt (dns* *2 + dew* *2) 
raddust = dustM* radmax 
hgtdust = dustM*dusthgt 

If (rad .It. 2 . 0* raddust .and. HGT .It. 2.0*hgtdust) 

& sizefact = 0.25* (1.0 + cos ( 90 . *DTR* rad/ raddust ) ) * 

Sc (1.0 + cos (90 . *DTR*HGT/hgtdust) ) 

Endif 

Computes the final dust storm factor values dustM and dustA, including the effects 
of time from start of the storm, the size factor based on position within the storm, 
and the storm intensity (0-3) . 


DSTF 34 
DSTF 35 
DSTF 36 
DSTF 37 
DSTF 38 
DSTF 39 
DSTF 40 
DSTF 41 
DSTF 42 
DSTF 43 
DSTF 44 
DSTF 45 


DSTF 27 
DSTF 28 
DSTF 29 
DSTF 30 
DSTF 31 
DSTF 32 
DSTF 33 


DSTF 20 
DSTF 21 
DSTF 22 
DSTF 23 
DSTF 24 
DSTF 25 
DSTF 26 


dustM = sizefact*dustM*intens/3 . 
dustA = sizefact*dustA*intens/3 . 


DSTF 46 
DSTF 47 
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Subrout ins : Da t a s t ep 


Code : DSTP 


Description: Driver routine for calling the atmospheric models and 
evaluating the winds and the density perturbations. 

Called By: MarsGRAM (Main) , Marsgrmb (Main) 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

COSPARNH, DATACOM, FILENAME 

Input Variables (not passed through commons): 

Description 


current height 
current latitude 
current West longitude 
current time 

Julian date for initial position 
height displacement between successive positions 
latitude displacement between successive positions 
longitude displacement between successive positions 
time displacement between successive positions 
flag for end of file on trajectory input 
position counter index 

previous density perturbation value for perturbation 
model 


Name Type 


CHGT .... 

REAL *4 

CLAT .... 

REAL *4 

CLON. . . . 

REALM 

CSEC .... 

REAL *4 

DATEO . . . 

REAL* 4 

DELHGT . . . 

REAL* 4 

DELLAT . . . 

REAL* 4 

DELLON . . . 

REAL *4 

DELTIME . . 

REAL *4 

EOF .... 

INTEGER* 4 

I 

INTEGER* 4 

RHO .... 

REAL *4 


Output Variables (not passed through commons): 


Name 

Type 

DENS .... 

REAL *4 

DENSHI . . . 

REAL *4 

DENSLO . . . 

REAL* 4 

DENSP . . . 

REAL *4 

EWWIND . . . 

REAL* 4 

NSWIND . . . 

REAL *4 

PRES .... 

REAL *4 

TEMP .... 

REAL* 4 


Description 


average atmospheric density 

high density (approx. +1 sigma) 

low density (approx. -1 sigma) 

density perturbation (% from mean value) 

Eastward wind component 

Northward wind component 

atmospheric pressure 

atmospheric temperature 


Variables (not passed through c ommo ns ) : 

Name Type Description 


ABSWLAT 
ALOGDENS 
ALS . . 

AMPl . . 

AMP 2 . . 

AREAHGT 
BRUNTF . 
BVF . . 

CD. . . 

CHGTS . 
CLATM . 
CL ATP . , 


REALM 
REAL *4 
REAL *4 
REALM 
REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL* 4 
REAL *4 
REAL* 4 
REAL *4 


absolute value of latitude, for wind calculations 

log (base 10) density 

areocentric longitude of the sun (Ls) 

diurnal amplitude in Zurek wave model 

semi-diurnal amplitude in Zurek wave model 

area-averaged terrain height 

Brunt-Vaisala frequency 

dummy argument in ATMOS 2 for Brunt-Vaisala frequency 
geostrophic drag coefficient 
height above local terrain 

negative displaced latitude for wind calculations 
positive displaced latitude for wind calculations 
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CLONM . . . 

REAL* 4 

CLONP . . . 

REAL*4 

CORIOL . . . 

REAL*4 

CORREL . . . 

REAL*4 

COSFAC . . . 

REAL* 4 

COSWLAT . . 

REAL*4 

DAMPING . . 

REAL*4 

DATE .... 

REAL* 8 

DCOS .... 

REAL*4 

DELEW . . . 

REAL*4 

DELNS . . . 

REAL*4 


DELZ .... 

REAL* 4 

DENO .... 

REAL *4 

DENOM . . . 

REAL* 4 

DENSURF . . 

REAL* 4 

DENSWA . . . 

REAL *4 

DEVAV . . . 

REAL* 4 

DEVHI . . . 

REAL *4 

DEVLO . . . 

REAL*4 

DHGT .... 

REAL*4 

DLATM . . . 

REAL* 4 

DLATP . . . 

REAL*4 

DLONM . . . 

REAL*4 

DLONP . . . 

REAL* 4 

DMINUS . . . 

REAL*4 

DPLUS . . . 

REAL *4 

DUSTA . . . 

REAL*4 

DUSTM . . . 

REAL*4 

FACTHI . . . 

REAL *4 

FACTLO . . . 

REAL* 4 

FACTOR . . . 

REAL *4 

FHI .... 

REAL*4 

FLO .... 

REAL *4 

FUG .... 

REAL *4 

FUGP .... 

REAL *4 

FVG .... 

REAL*4 

FVGP .... 

REAL* 4 

GZ 

REAL *4 

HLATM . . . 

REAL *4 

HLATP . . . 

REAL* 4 

HLONM . . . 

REAL*4 

HLONP . . . 

REAL*4 

HLS .... 

REAL *4 

HSCALE . . . 

REAL*4 

IFAULT . . . 

INTEGER* 4 

L 

INTEGER* 4 

MARSAU . . . 

REAL*4 

OHGT .... 

REAL*4 

OHGTS . . . 

REAL*4 

PCOS .... 

REAL*4 

PERT .... 

REAL*4 

PERTMAX . . 

REAL* 4 

PERTMIN . . 

REAL*4 

PI180 . . . 

REAL*4 

PLATM . . . 

REAL*4 

PLATP . . . 

REAL*4 

PLONM . . . 

REAL*4 

PLONP . . . 

REAL* 4 

RLATM . . . 

REAL* 4 

RLATP . . . 

REAL* 4 

RLONM . . . 

REAL*4 

RLONP . . . 

REAL*4 

RSC .... 

REAL*4 


negative displaced longitude for wind calculations 
positive displaced longitude for wind calculations 
Coriolis parameter for wind calculations 
one-step correlation value for density perturbations 
cosine factor used in wind calculations 
cosine of absolute value of latitude 

viscous damping factor in mountain wave perturbations 

current Julian date 

COSPAR model atmospheric density 

E-W displacement between previous and current 

position (DSTP 77); also E-W distance (= 5° 
longitude), for wind calculation (DSTP129) 

N-S displacement between previous and current 
position (DSTP 76); also N-S distance (= 5° latitude), 
for wind calculation (DSTP120) 
height displacement, previous to current position 
dummy surface density argument for ATM0S2 subroutine 
denominator term in wind calculations 
atmospheric density at surface 
density wave perturbation amplitude, % 
deviation of average density from COSPAR value, % 
deviation of high (+1 sigma) density from COSPAR, % 
deviation of low (-1 sigma) density from COSPAR, % 
terrain height factor for mountain wave magnitudes 
dummy density argument used in calling ATMOS2 

dummy density argument used in calling ATM0S2 

dummy density argument used in calling ATMOS 2 

dummy density argument used in calling ATMOS 2 

difference between average and low density 
difference between high and average density 
relative dust storm magnitude for diurnal variations 
relative dust storm magnitude for daily mean values 
ratio of high (+1 sigma) to average density 
ratio of average to low (-1 sigma) density 
boundary layer wind factor 

dummy high density factor for call to ATM0S2 

dummy low density factor for call to ATM0S2 

intermediate E-W component in wind calculations 

intermediate E-W component in wind calculations 

intermediate N-S component in wind calculations 

intermediate N-S component in wind calculations 

acceleration of gravity at current height 

dummy scale height argument for call to ATMOS 2 

dummy scale height argument for call to ATM0S2 

dummy scale height argument for call to ATM0S2 

dummy scale height argument for call to ATM0S2 

horizontal scale for density perturbations 

pressure scale height 

error flag for the PPND function 

error flag for the RANDOM function 

Sun-Mars orbital radius at current time 

height above reference ellipsoid (for output purposes) 

height above local terrain (for output purposes) 

COSPAR model atmospheric pressure 
mountain wave perturbation magnitude 

maximum perturbation magnitude allowed by stability 
minimum allowed perturbation magnitude 
factor to convert degrees to radians 
pressure south of position, for pressure gradients 
pressure north of position, for pressure gradients 
pressure east of position, for pressure gradients 
pressure west of position, for pressure gradients 
dummy radius from Mars center, used in ATM0S2 call 

dummy radius from Mars center, used in ATM0S2 call 

dummy radius from Mars center, used in ATMOS 2 call 

dummy radius from Mars center, used in ATMOS 2 call 

radius from center of Mars to current position 


A-24 



SIGD. . 
SINWLAT 
SUNLAT. 
SUNLON . 
TAVG . . 
TAVM. . 
TAVP . . 

TCOS . . 

TFAC . . 

TFACTOR 
THGT . . 

TLATM . 
TLATP . 
TLOCAL . 
TLONM . 
TLONP . 
TMAM. . 
TMAP . . 

TMAX. . 

TMIM. . 

TMIN . . 

TMIP . . 

VAR . . 

VARX. . 

VARY. . 

vise. . 

VISCFAC 
VLL . . 

VLS . . 

WAVE. . 
WIND . . 

WLAT. . 

2 . . . 
20 . . . 
21 . . . 
22 . . . 
ZF. . . 

Z FACTOR 
ZT. . . 


REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL * 4 
REAL *4 
REALM 
REAL *4 
REAL* 4 
REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REALM 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REALM 
REAL *4 
REAL *4 
REAL *4 
REALM 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 


Methodology: 


standard deviation for perturbations {% of mean) 

sin of absolute value of latitude 

latitude of sub-solar position 

West longitude of sub-solar position 

daily average surface temperature 

dummy average surface temperature in ATMOS 2 call 

dummy average surface temperature in ATMOS 2 call 

COSPAR model atmospheric temperature 

dummy density perturbation factor for call to ATMOS 2 

density perturbation factor at base of thermosphere 

local terrain height 

dummy temperature value for call to ATMOS 2 

dummy temperature value for call to ATMOS 2 

local time in "Martian hours" (l/24th Sols) 

dummy temperature value for call to ATMOS 2 

dummy temperature value for call to ATMOS 2 

dummy maximum surface temperature in ATMOS 2 call 

dummy maximum surface temperature in ATMOS 2 call 

daily maximum surface temperature 

dummy minimum surface temperature in ATMOS 2 call 

daily minimum surface temperature 

dummy minimum surface temperature in ATMOS 2 call 

output height variable for the OUTPUT file 

x (1st) variable for plotable output files 

y (2nd) variable for plotable output files 

coefficient of molecular viscosity 

viscous modification factor for winds 

vertical scale for wind viscosity factor 

ver tical scale density for perturbations 

relative density perturbation from Zurek wave model 

wind for use in perturbation magnitude 

latitude used for wind calculations 

height above ground (m) , for boundary layer winds 

surface roughness (m) 

random variable, used in perturbation calculations 

random variable, used in perturbation calculations 

height of base of thermosphere 

height factor in vertical interpolation 

dummy thermosphere base height, used in ATMOS 2 call 


Datastep is the principal driver subroutine for stepping through the 
positions and evaluating the atmospheric variables. It starts by reading 
the next position from the trajectory input file or by automatically 
generating the next position from the increments in latitude, longitude 
altitude and time. 

Terrain height is obtained by calling the Terrain subroutine (DSTP 
34), orbit parameters are evaluated by calling the ORBIT subroutine (DSTP 
36), dust storm perturbation factors (dustM and dustA) are evaluated by 
calling subroutine Dustfact (DSTP 38) and the other mean atmospheric 
parameters are evaluated by calling ATMOS 2 (DSTP 40) . High (+1 sigma) and 
low (-1 sigma) density values are computed from the factors FACTHI and 
FACTLO (DSTP 43 - DSTP 57). A normally-distributed random variate is found 
using RANDOM and PPND (DSTP 58 - DSTP 61). The vertical scale for the wind 
viscosity factor is set to the pressure scale height (DSTP 65). The 
vertical and horizontal scales for the density perturbations are found from 


VLS =8.0 


DSTP 69 


HLS = 30. + . 01875*CHGT**2 


DSTP 72 
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IF {HLS . GT . 600 . ) HLS = 600. 


DSTP 73 


The components of the relative displacement between previous and current positions 
are computed 


DELNS = DTR*RSC* (DELLAT) /HLS 

DELEW = -DTR*RSC*COS (DTR*CLAT) *DELLON/HLS 

DELZ = DELHGT/VLS 


DSTP 76 
DSTP 77 
DSTP 78 


and used to evaluate the correlation between the previous and current position 


CORREL = Exp (- (Abs{ DELNS) + Abs (DELEW) + Abs(DELZ))) 


DSTP 81 


Next, the wind components are computed from the horizontal pressure gradients 
(a Mars version of the geostrophic wind relations, referred to as the areostrophic 
winds; see discussion on computation of winds in the section "WIND PROFILES IN MARS- 
GRAM" , of the Release 1 Technical Report, 1989; see Appendix B) . The pressure 
gradients are computed by finite differences of the pressure across +/- 2.5 
latitude and longitude displacements from the current position, by calling ATMOS 2 
(longitude + 2.5° at DSTP 86; longitude - 2.5° at DSTP 91; latitude + 2.5° at DSTP104 
and latitude -2.5° at DSTP 113). The area-average terrain height, for use in the 
mountain-wave perturbation model, is evaluated with the Terrain function (DSTP 95 
and DSTP117) . For computing the pressure gradients, the distances corresponding to 
the 5 degree latitude or longitude displacements are computed at DSTP120 and 
DSTP129. The areostrophic wind components, computed from the pressure gradients are 
evaluated by 


FUG = - { PLATP- PLATM) / (DELNS* DENS) 
FVG = ( PLONP-PLONM) / (DELEW 1 *' DENS) 


DSTP137 

DSTP138 


[Release 1 Technical Report, equation (1); see Appendix B] . 

For modifications by molecular viscosity (relevant in the thermospheric 
region) , viscosity factors are computed 


VISC = BETA*TEMP* *1.5/ (TEMP + SVAL) 
VISCFAC = VISC/ (1.0E6*DENS*VLL**2) 


DSTP133 

DSTP134 


For absolute latitudes less than 7.5°, a special low-latitude wind model 
[Release 1 Technical Report, equations (2) and (3); see Appendix B] is used 


FUGP = -4000 . *RSC*SINWLAT* ( PLATP + PLATM - 2.*PRES)/ 
& { DENS * COSWLAT * DELNS * DELNS ) 

FVGP =0.0 


DSTP141 

DSTP142 

DSTP143 


For absolute latitudes between 7.5° and 15°, winds are interpolated between the low 
latitude and regular values 


FUG = FUGP + (ABSWLAT-7 . 5) * (FUG-FUGP) /7 . 5 
FVG = FVGP + (ABSWLAT-7. 5) * (FVG-FVGP) /7 . 5 


DSTP145 

DSTP146 
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For absolute latitudes greater than 75°, a special high-latitude relation [Release 1 
Technical Report, equation (4); see Appendix B] is used 


Cosfac = COSWLAT/Cos (pil80*75 . ) DSTP154 

FUG = FUG*COSWLAT*Cosfac DSTP155 

FVG = FVG*COSWLAT*Cosfac DSTP156 

Viscous-corrected wind components are computed by [Release 1 Technical 
Report, equations (6a) and (6b); see Appendix B] 

DENOM = CORIOL* * 2 + VISCFAC**2 DSTP158 

EWWIND = ( CORIOL*FUG - VISCFAC* FVG) /DENOM DSTP160 

NSWIND = (CORIOL*FVG + VISCFAC* FUG) /DENOM DSTP161 

altitudes within 1 km of the surface, boundary— layer modifications are appl: 
[Release 1 Technical Report, equations (8) and (9); see Appendix BJ 

factor = ( Sqrt (CD) /0.4)*Alog(z/z0) DSTP172 

if (factor . gt. 1.0) factor = 1.0 DSTP173 

EWWIND = factor *EWWIND DSTP174 

NSWIND = factor * NSWIND DSTP175 

where z is the height above the surface (m) , the drag coefficient is taken to be 
= 0.0015 and the surface roughness is assumed to be zO = 0.03 m. 

For the mountain wave perturbation model [Release 1 Technical Report, 
equations (10) and (15); see Appendix B] , the steps are: 

(1) Compute the terrain height adjustment factor dhgt (local height above area- 
averaged terrain height, m) by 

dhgt = 1000.*(thgt - areahgt) DSTP184 

dhgt — 2 . 0 *dhgt DSTP187 

If ( dhg t .It. 10.) dhgt = 10. DSTP1 8 8 

If (dhgt .gt. 2000.)dhgt = 2000. DSTP189 

(2) Compute the wind speed 

Wind = Sqrt (NSWIND**2 + EWWIND**2) DSTP191 

If (Wind .It. 1.0) Wind = 1.0 DSTP192 

(3) Compute the viscous damping factor 

Damping = 1 . 0E6 *HLS*VISC* (Brunt f * * 3 ) *Hscale* DSTP194 

& (1. - (DENS/densurf ) ) / (12 . 5 663 7 *DENS*Wind* *4 ) DSTP195 

If (Damping . gt . 50.) Damping = 50. DSTP196 

(4) Compute the perturbation magnitude and insure that it does not exceed the 
maximum magnitude allowed by stability and that it is at least as large as a 
prescribed minimum value (dependent on height) 
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DSTP198 


pert = FACTHI*dhgt*Exp (-Damping) 

pertmax = (500 . *VLS*Bruntf**2/gz) * (1 . + 0 . 5*VLS/Hscale) 

If (pertmax .gt. 0 . 5*VLS/Hscale) pertmax = 0 . 5*VLS/Hscale 

If (pert .gt. pertmax) pert = pertmax 

pertmin = 0.01 + 0.001*CHGT 

If (pertmin .gt. 0.05) pertmin =0.05 

If (pert .It. pertmin) pert = pertmin 


DSTP199C 

DSTP199d 

DSTP200 

DSTP201 

DSTP202 

DSTP203 


(5) Between 75 km and the base of the thermosphere, interpolate between the mountain 
wave model and the magnitude derived from the Stewart thermosphere model at the base 
of the thermosphere 


zf actor = (CHGT - 75.0) /(ZF - 75.0) 
pert = pert + zfactor* ( tfactor - pert) 


DSTP207 

DSTP208 


(6) Values of DENSHI and DENSLO are computed, based on the mountain wave 
perturbation magnitudes 


585 DENSHI = DENSM1. + pert) 
DENSLO = DENS/{1. + pert) 
DPLUS = DENSHI - DENS 
DMINUS = DENS - DENSLO 


DSTP212 

DSTP213 

DSTP214 

DSTP215 


The local time, in Martian hours = l/24th Sols is computed from the longitude 
of the current position 


T LOCAL = 12. + (SUNLON - CLON)/15. 

IF ( TLOCAL .LT. 0.)TLOCAL = TLOCAL + 24. 
IF (TLOCAL .GT. 24.) TLOCAL = TLOCAL - 24. 


DSTP218 

DSTP219 

DSTP220 


Height above the reference ellipsoid or above the local surface terrain 
height are evaluated and saved for output purposes 


OHGT = CHGT 
OHGTS = CHGTS 
IF (OHGT .LE. -5.) THEN 
OHGT = thgt 
OHGTS = 0 . 

END IF 


DSTP223 

DSTP224 

DSTP225 

DSTP226 

DSTP227 

DSTP228 


The Zurek wave perturbation model 
Observations and Model Profiles for Mars 
Pitts et al., eds . , is used to add tidal 


(Section 4.3 of "The Mars Atmosphere: 
Missions", Report No. JSC-24455, David E. 
wave perturbations 


If (modpert .ne. lJCall Wavepert (OHGT, CLAT, TLOCAL, Dus tM, Dus tA, 
& wave , ampl , amp2 ) 


DSTP231 

DSTP232 


The Zurek tidal wave model and the mountain-wave model perturbation are added, with 
the mountain wave values treated as having the previously-computed correlation 
between the prior and current positions 


RHO = CORREL*RHO + SQRT (1.0 - CORREL** 2 ) * Zl 

IF (RHO. LT . 0 . 0) DENSP = DENSM1. + wave) + RHO*DMINUS 


DSTP235 

DSTP236 
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IF {RHO . GE . 0 . 0 ) DENSP = DENSM1. + wave) + RHO*DPLUS 
If (DENSP .It. 0.0) DENSP = 0.05*DENS 


DSTP237 

DSTP238 


The standard deviation and DENSHI and DENSLO values are adjusted to reflect 
the Zurek wave perturbation values 


SIGD = 50. * (DENSHI -DENSLO) /DENS 


DSTP241 


DENSHI = DENSHI + (ampl + amp2)*DENS 
DENSLO = DENSLO - (ampl + amp2)*DENS 


DSTP243 

DSTP244 


DENSP = 100.* (DENSP - DENS) /DENS 


DSTP246 


Deviations from the COSPAR reference atmosphere are computed (DSTP254b to 
DSTP254k) . Descriptively formatted data are written to the LIST file ( DSTP248 and 
DSTP255) . Output to the various plotable-output (graphics-input) files (units 21- 
34) are written, according to the options selected (NVARX, NVARY , and logscale; 
DSTP264b through DSTP349) . y 
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Subroutine: DZDUST 


Code : DZDS 


Description: calculates the perturbation to ZF, the height of the base of 
the thermosphere, due to dust storms. Modified from the original Stewart 
subroutine, which assumed storms starting at Ls = 205 and Ls = 275 each 
year. Now uses the same starting Ls value as for lower atmosphere (if a 
dust storm case is selected) . 


Called By: ATMOS 2 , ST EWART 2 

Common Blocks used (See Table 2-3 for a list of variables in c 

None 

Input Variables (not passed through commons) : 

Name Type Description 


blocks) : 


INTENS. 

. . REAL *4 

ALSO. . 

. . REAL*4 

LSUN. . 

. . REAL* 4 

Output 

Variables (not 

Name 

Type 

DZ. . . 

. . REAL *4 

Local 

Variables (not 

Name 

Type 

QA. . 

. . . REALM 

DZA . 

. . . REAL *4 

DLS . 

. . . REAL* 4 

DLSA. 

. . . REALM 


dust storm intensity level (0.0 - 3.0) 
starting Ls value for storm 
areocentric longitude of sun (Ls) 


Description 


increase in height of base of thermosphere, due to 
dust storm, Stewart model method 


Description 


exponential decay factor for dust storm effect 
intermediate value for DZ calculation 
Ls duration for dust storm effect (25 degrees of Ls) 
difference in Ls value from start of storm (Ls - LsO) 


Methodology : 

This subroutine has been modified from the original subroutine of 
Stewart (which still appears in the code, but commented out) . Instead of 
assuming two storms per Mars year (one starting at Ls = 205, another 
starting at Ls = 275), DZDUST now uses the same Ls value for the beginning 
of the storm as is used for the lower altitude regions (based on user input 
value for LsO) . 

The Ls difference, Ls-LsO, is put in the range 0-25, by the relation 


DLSA = AMOD ( ( ( LSUN - alsO) + 720.0), 360.0) 


The exponential time-decay factor is computed by 


QA = EXP ( -DLSA / DLS) 


DZDS 28 


A-30 


and the final displacement for the base of the thermosphere (km) is computed by 


DZA = 5 . 0*INTENS 

DZ = DZA*QA* (1.0 - QA* *4 ) *1.869 


DZDS 29 
DZDS 30 
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Sub rout in* : GScalc 


Cod* : ESCL 


Description: Computes short-term and long-term standard deviations from 
nominal densities in Stewart thermosphere model. 

C*ll*d By: ATMOS 2, STEWART 2 

Coma on Blocks Used (S*« Tabl* 2-3 for » list of variables in common blocks): 

None 

Input Variables (not p****d through commons) : 

Name Type Description 


SIGMA. . . REAL*4 long-term thermospheric variability factor 

STDL . . . REALM short-term thermospheric variability factor 

Output Variable* (not passed through coonons) : 

Name Type Description 


gg .... REAL*4 array for thermospheric variability factors 


T rfwi Variable* (not passed through cossoons) : 

Name Type Description 


j INTEGERM index for the ES array 

EPS. . . . REAL*4 intermediate array for ES calculations 

gX(j . . . . REAL*4 array of coefficients for ES calculations 


Methodology: 

The short-term thermospheric variability factor, stdl, is read in the 
main (MarsGRAM) or SETUP routines. The long-term variability factor, 

SIGMA, is set by the calling value in the STEWART2 argument list (e.g. 0.0 
at ATM2143 , 1.0 at ATM2146 and ATM2164, or -1.0 at ATM2166) . 

Even-number indexed ES values (0, 2, 4, 6, 8, 10) are for long-term 
variability; odd-number indexed ES values (1, 3, 5, 7, 9, 11) are for 
short-term variability. ES(0) and ES(1) are for the dependence on solar 
activity ( FBAR at STW2 25). ES(2) and ES<3) are for the exospheric 

temperature (TINF at STW2 46). ES(4) through ES(7) are for the atomic 

oxygen factors (FO at THRM 38 and AO at THRM 36). ES(8) and ES(9) are for 

the base-height of the thermosphere (ZF at ATM2137 and STW2 38) and the 
temperature at the base of the thermosphere (TF at STW2 47). ES(10) and 

gg^Xl) are for the dust storm effect on the thermosphere (DUST at ATM2135 

and STW2 37) . , 

The factors ES(1), ES(6) and ES(ll) are set to zero m the Stewart 
thermosphere model, and have therefore been left out of the equations that 
they would otherwise appear in (i.e. they do not appear explicitly in lines 
STW2 25 [ES ( 1 ) ] , THRM 36 [ES(6)] or in ATM2135 or STW2 37 [ ES ( 11 ) ] . 

Even-number indexed values of EPS are set to stdl; odd-number indexed 
values of EPS are set to SIGMA. The SIG array elements are set to assigned 
values and ES is computed as the product of EPS and SIG. 
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Subroutine: gamma 


Code: GAMA 


Description: Computes lapse rates (-dT/dz) over the significant level 
areopotential height intervals, for a given latitude and Ls value. 

Called By: ATMOS 2 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

None 


Input 

Variables (not 

passed through commons) : 

Name 

Type 

Description 

ALAT. 

. . . REAL *4 

current latitude 

ALS . 

. . . REAL* 4 

areocentric longitude of the sun (Ls) 

Output Variables (not passed through commons) : 

Name 

Type 

Description 

GAM . 

. . . REAL* 4 

lapse rate array (negative of temperature gradient) 

Local 

Variables (not 

passed through commons) : 

Name 

Type 

Description 


A . . 
Al . . 
A2 . . 
A3 . . 
B . . 
B1 . . 
B2 . . 
B3 . . 
C . . 
Cl. . 
C2 . . 
C3 . . 
CALS. 
F180. 
F240 . 
I . . 
PI120 
PI180 
PI90 . 
SALS. 


REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL*4 
INTEGER* 4 
REAL *4 
REAL* 4 
REAL* 4 
REAL* 4 


coefficient for part of gamma depending only on latitude 

gam { 1 ) (0-5 km) coefficients for A term 

gam(2) (5-15 km) coefficients for A term 

gam(3) (15-30 km) coefficients for A term 

coefficient depending on latitude and sine of Ls 

gam(l) (0-5 km) coefficients for B term 

gam(2) (5-15 km) coefficients for B term 

gam(3) (15-30 km) coefficients for B term 

coefficient depending on latitude and cosine of Ls 

gam(l) (0-5 km) coefficients for C term 

gam(2) (5-15 km) coefficients for C term 

gam(3) (15-30 km) coefficients for C term 

cosine of Ls angle 

array of Fourier terms for tc/ 90 times latitude 

array of Fourier terms for 7C/120 times latitude 

index value 

rt /120 

71/180 

71/90 

sine of Ls angle 


Methodology : 

Uses Fourier coefficients for (sine and cosine) terms in latitude and 
Ls angle. The temperature lapse rate (-dT/dz) array, gam (in units of 
K/km) , has 5 elements, one for each height interval between the significant 
levels: gam(l) applies to heights 0-5 km, gam(2) is for 5-15 km, gam(3) is 
for 15-30 km, gam(4) is for 30-50 km, and gam{5) applies to 50-75 km. The 
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upper two intervals are assumed to have fixed lapse rates: gam(4) = 0.9 
K/kro and gam(5) =0.4 K/km. For the three lower levels, the gam values are 
computed by an equation of the form 


gam = A + B sin(Ls) + C cos(Ls) 


(GAMA 58, GAMA 68 and GAMA 78) . 

The A coefficient term is computed as a Fourier series, in terms of 
sines and cosines of multiples of n/90 times the latitude (GAMA 54, GAMA 64, and 
GAMA 74). The B and C coefficient terms are computed as Fourier series, in 
terms of sines and cosines of multiples of it/120 times the latitude (GAMA 55, 
GAMA 57, GAMA 65, GAMA 67, GAMA 75 and GAMA 77). Coefficients (Al, A2 , A3, B1 , B2 
g3 ' Cl, C2 , and C3 arrays) are stored in data statements, indexed 0 through 8, fo^ 1 
the constant coefficient and the 8 Fourier terms in the respective latitude series. 
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Subroutine: orbit 


Code : ORBT 


Description: Computes orbital radius of Mars from Sun, areocentric 
longitude of the sun (Ls) and the latitude and longitude of the sub-solar 
point on the Mars surface for a given time (expressed as a Julian date) . 

Called By: MarsGRAM (Main) , Datastep, SETUP 

Common Blocks Used (See Table 2-3 for a li.t of variables in common blocks) : 

None 


Input Variables (not passed through commons) : 


Name Type 


Description 


XDATE . 

. . REAL* 8 

Julian date for current time 


Output 

Variables (not passed through commons): 


Name 

Type 

Description 


LATSUN . 
LONSUN . 
LSUBS . 
RADIUS . 

. . REAL *4 
. . REAL* 4 
. - REAL* 4 
. . REAL *4 

latitude of sub-solar point on Mars on Mars surface 
West longitude of sub-solar point on Mars surface 
areocentric longitude of the sun (like right ascension) 
magnitude of vector from Sun to orbital position of Mars 

Local Variables (not 

passed through commons) : 


Name 

Type 

Description 



AO thru 

A6 

REAL *8 

BO thru 

B6 

REAL *8 

CO thru 

C6 

REAL* 8 

DO thru 

D6 

REAL *8 

DATE . . 


REAL *8 

DS . 


REAL *8 

LON . . 


REAL *8 

LONO . . 


REAL *8 

LS. . . 


REAL *8 

LSO . . 


REAL* 8 

PERI . . 


REAL *8 

PER2 . . 


REAL* 8 

PERLON . 


REAL *8 

PERLS . 


REAL *8 

PI160 . 


REAL *8 

RAD . . 


REAL* 8 

TIME1 . 


REAL* 8 

TIME2 . 


REAL *8 


TWOPI 
XLON . 
XLS . 


REAL* 8 
REAL *8 
REAL* 8 


Fourier coefficients for RADIUS calculation 
Fourier coefficients for LATSUN calculation 
Fourier coefficients for LSUBS calculation 
Fourier coefficients for LONSUN calculation 
modified Julian date (xdate - 2442779. + 0.5) 
intermediate value of LATSUN used for calculation 
intermediate value of LONSUN used for calculation 
initial value of LONSUN 

intermediate value of LSUBS used for calculation 
initial value of LSUBS 

687 day period (used for LATSUN, LSUBS and RADIUS) 
696 day period (used for LONSUN) 

coefficient for linear increase of LONSUN with date 
coefficient for linear increase of LSUBS with date 
n / 180 

intermediate value of RADIUS used for calculation 
current time variable with period = PERI 
current time variable with period = PER2 
2 n 

part of LONSUN increasing linearly with date 
part of LSUBS increasing linearly with date 


Me t hodol ogy : 


r*momT USeS Pre com P ut:ed coefficients in Fourier time series to compute 
TSUN, LONSUN, LSUBS and RADIUS from the modified Julian date 
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DATE = XDATE - 2.442779d6 + 0 . 5d0 


ORBT 18 


The parts of LSUBS and LONSUN that increase linearly with the date are 
computed from 


xls = IsO + 3 . 60d2*date/perls 

xlon = lonO + 3.60d2Mdate - 2 . 922d3 ) /perlon 


ORBT 19 
ORBT 20 


The time variables for period 1 (687 day period) and period 2 (696 day period) are 
evaluated by 


TIMEl=DATE/PERl 
TIME2 =DATE/ PER2 


ORBT 21 
ORBT 22 


LSUBS (intermediate value LS) , in degrees, is computed from 


LS=C0+C1*DSIN (TIME1 ) +C2*DCOS (TIMEl ) +C3 *DSIN (2 . 0d0*TIMEl) 

& +C4*DCOS ( 2 . OdO*TIMEl ) +C5*DSIN ( 3 . OdO*TIMEl) +C6*DCOS ( 3 . OdO*TIMEl ) 
& +XLS 


ORBT 23 
ORBT 24 
ORBT 25 


LATSUN (intermediate value DS) , in degrees, is computed from 

DS=B0+B1*DS IN (TIMEl ) +B2*DCOS (TIMEl ) +B3*DSIN(2 . 0d0*TIMEl) ORBT 26 

& +B4*DCOS (2 . 0d0*TIMEl ) +B5*DSIN ( 3 . 0d0*TIMEl) +B6*DCOS (3 . 0d0*TIMEl ) ORBT 27 


LONSUN (intermediate value LON), in degrees, is computed from 


LON=DO+D1*DSIN (TIME2 ) +D2*DCCS (TIME2 ) +D3*DSIN (2 . 0d0*TIME2 ) ORBT 28 

& +d4*DCOS(2.0d0*time2)+d5*DSIN(3.0d0*time2)+d6*DCOS(3.0d0*time2) ORBT 29 

n ORBT 3 0 

& +xlon 

RADIUS (intermediate value RAD), in astronomical units (AU) , is computed from 


rad = aO + al*DSIN(timel) + a2*DCOS (timel) 

& + a3*DSIN(2 . 0d0*timel) + a4*DCOS (2 . 0d0*timel) 
& + a5*DSIN (3 . 0d0*timel) + a6*DCOS (3 . 0d0*timel) 


ORBT 35 
ORBT 36 
ORBT 37 


The Fourier coefficients (A0-A6, B0-B6, C0-C6 and D0-D6) are loaded in via 
data statements. Both LSUBS and LONSUN are put into the range 0 to 360 
degrees . 
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Function: 


phasint 


Cod® : PHSN 


Description: Linear interpolation on phase angle, by vector components, for 
the phases of the wave perturbation field. 

Called By: Wavepert 

Common Blocks Used (See Table 2-3 for a list of variables in comd blocks): 

None 

Input Variables (not passed through commons) : 

Description 


linear interpolation value (0-1) 

frequency for the phase angle interpolation 

phase angle at first position (location where DX=0) 

phase angle at second position (location where DX=1) 


Name Type 


DX REALM 

OMEGA . . . REAL *4 

PHI .... REAL* 4 
PH2 .... REALM 


Output Variables (not passed through commons) : 


Name 

Type 

Description 




PHASINT . . REAL *4 

interpolated phase angle 

at time 

corresponding input 

Local 

Variables (not 

passed through con 

naons) : 



Name 

Type 

Description 




DPH . 

. . . REAL* 4 

phi2 - phil 




X . . 

. . . REAL* 4 

interpolated "x" 

vector component 

of phase 

Xl . . 

. . . REAL* 4 

"x" component of 

phase at 

time 1 

[sin (omega*phil ) ] 

X2 . . 

. . . REAL *4 

"x" component of phase at 

time 2 

[sin(omega*phi2 ) ] 

Y . . 

. . . REAL* 4 

interpolated "y" 

vector component 

of phase 

Yl . . 

. . . REAL* 4 

"y" component of phase at 

time 1 

[cos (omega*phil) ) 

Y2 . . 

. . . REAL* 4 

"y" component of phase at 

time 2 

[cos (omega*phi2 ) ] 


Methodology : 

This routine is used to interpolate wave phase angle as a function of 
latitude or height, in the Zurek wave model (subroutine Wavepert) . It 
treats phase angle as a circular variable, having "x" and "y" Cartesian 
vector components . The input parameter omega is the circular frequency for 
the latitude or height variation. Input values of phase angle (phil and 
phi2) corresponding to x and y component values 


xl = sin (omega*phl) 
x2 = sin (omega*ph2 ) 
yl = cos (omega*phl) 
y2 = cos (omega*ph2 ) 


The interpolated values of the vector 


x = xl + (x2 - xl ) *dx 
y = yl + (y2 - yl)*dx 


PHSN 4 
PHSN 5 
PHSN 6 
PHSN 7 


components are 


PHSN 8 
PHSN 9 
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The interpolated value of the phase angle is reconstructed from the vector 
components by 


phasint = atan2 (x, y) /omega 


PHSN 10 


The routine is also used to interpolate linearly on 24 hour local time 
variables , by the relations 


dph = ph2 - phi 

If (abs(dph) .gt. 12.)dph = dph - sign(24. ,dph) 
phasint = phi + dph*dx 


PHSN 12 
PHSN 13 
PHSN 14 
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Function: polecap 


Coda: POLC 


Description: Computes radius (latitude degrees) of the edge of the polar 
cap from the pole, for a given areocentric longitude of sun (Ls) for 
either north polar cap or south polar cap. 

Called By: Alb, Tsurface 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

None 


Input 

Variables (not 

passed through commons) : 

Name 

Type 

Description 

ALAT. . 
ALS . . 

. . REAL *4 

. . REAL *4 

latitude (+ for north pole, - for south pole) 
areocentric longitude of the sun (Ls) 

Output 

Variables (not passed through commons) : 

Name 

Type 

Description 

POLECAP 

. - REAL *4 

polar cap radius (degrees of latitude) 

Local Variables (not 

passed through commons) : 

Name 

Type 

Description 

ALSO. . 
PI180 . 

. . REAL * 4 

. . REAL *4 

polar cap phase in Ls angle 
71 / 180 


Methodology : 

Uses a sine wave in Ls to approximate the latitude progression and 
regression of the polar caps. The radial distance (in degrees of latitude) 
from the pole to the cap edge is computed by 


polecap = 19. - 16 . *sin (pil80* (als-alsO) ) 
where alsO = 230 for the north polar cap and alsO 


POLC 10 


50 for the south polar cap. 
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Function: PPND 


Code: PPND 


Description: Produces a normally-distributed deviate, for use in the random 
perturbation model. 

Callad By: Datastep 

Com ton Block* U»»d (S«« Tabl* 2-3 Cor a li*t of variable* in common block*) s 


None 

Input Variables (not passed through commons) : 

Name Type Description 

p REALM input value, from uniformly distributed random number 

Output Variables (not passed through commons) : 

Name Type Description 

I FAULT . . . INTEGERM error flag 

PPND. . . . REAL* 4 normally (Gaussian) distributed random variable 


bocal Variables (not passed through commons): 


Name 


Type 

Description 





AO thru 

A3 

REAL *4 

coefficients 

used 

in 

computing 

PPND 

Bl thru 

B4 

REALM 

coefficients 

used 

in 

computing 

PPND 

CO thru 

C3 

REAL*4 

coefficients 

used 

in 

computing 

PPND 

D1 thru 

D2 

REAL *4 

coefficients 

used 

in 

computing 

PPND 

HALF. . 


REAL* 4 

1./2. 





ONE . . 


REAL *4 

1. 





Q . • . 


REAL *4 

p - 1./2. 





R . . . 


REAL*4 

the square of Q 




SPLIT . 


REAL* 4 

0.42 





ZERO . . 

* * 

REAL *4 

0 . 





Methodology 

. 







Uses the method of algorithm AS 111, Applied Statistics (1977), 

Volume 26, page 118. The empirical coefficients (AO through C3 ) are used 
to fit the normally-distributed value PPND corresponding to the lower tail 
area p of the normal error integral. The error flag ifault is set to 0 if 0 
< p < 1. Otherwise, ifault is set to 1 and PPND is set to 0. 
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Subroutine: Pressure 


Code: PRES 


Description: Computes an array of atmospheric pressures at the significant 
height levels, from the surface pressure, the terrain height, surface 
gravity, local radius of reference ellipsoid, and temperatures at the 
significant levels. Also computes the areographic altitudes, z(i), 
corresponding to the areopotential altitudes, gph(i) , for the significant 
levels. 

Called By: ATM0S2 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks): 

None 

Input Variables (not passed through commons) : 


Name 

Type 

Description 

GAM . . 

. . REALM 

array of lapse rates across the significant levels 

GZERO . 

. . REAL *4 

surface gravity 

PSURF . 

. . REAL* 4 

surface pressure 

RBAR . . 

. . REAL *4 

gas constant 

RREF . . 

. . REAL *4 

local radius of reference ellipsoid 

T . . . 

. . REAL* 4 

array of temperatures at the significant levels 

THGT . . 

. . REAL *4 

terrain height above reference ellipsoid 

Output 

Variables (not 

passed through commons) : 

Name 

Type 

Description 

P . . . 

. . REAL *4 

array of pressures at the significant levels 

z . . . 

. . REAL* 4 

array of areographic heights at the significant level: 

Local Variables (not passed through commons) : 

Name 

Type 

Description 

I . . . 

* . INTEGER* 

4 array index 

GOR . . 

- . REAL* 4 

gravity divide by gas constant 

DH. . . 

. . REAL* 4 

height increment array between significant levels 

GPH . . 

. . REAL *4 

areopotential height array of the significant levels 


Methodology : 

Evaluates the areographic altitudes, z(i), that correspond to the 
areopotential heights, gph(i), of the significant levels (0, 5, 15, 30, 50 
and 75 areopotential km). 


using 


z ( i ) “ ( 9Ph ( i ) + thgt ) / ( 1 . - (gph ( i ) +thgt) /rref ) PRES 22 

Evaluates the pressure array p(i) at the significant level altitudes, 


P(i) = P(i-l) *exp{-goR*alog(T{i-l) /T(i) ) /gam(i) ) 


PRES 25 
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if the lapse rate, gam(i), is not zero, or 

p(i) = p ( i— 1 ) *exp (-goR*dh ( i) /T ( i-1 ) ) 

if the lapse rate is zero. 


PRES 27 
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Subroutine: PRSEAS 


Code: PRSE 


Description: Computes the relative seasonal variation in atmospheric 

pressure on the reference ellipsoid, at a given latitude and areocentric 
longitude of the sun (Ls) . 

Called By: ATMOS 2 , Psurface, STEWART2 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

None 

Inpub Variablds (nob passed bhrousfh commons) : 

Description 

latitude of current position 
areocentric longitude of the sun (Ls) 

passed bhrough commons) : 

Description 

relative pressure variation with latitude and season at 
height of reference ellipsoid 


Name Type 

DAT .... REAL *4 

LSUN. . . . REAL* 4 

Oubpub Variables (nob 

Name Type 

PR REAL *4 


Local Variables (nob passed bhrough commons) : 


Name Type 


Al . . 
All . 
A12 . 
A2 . . 
A21 . 
A22 . 
LAT2 . 
PHI 1 . 
PHI11 
PHI 12 
PHI 2 . 
PHI21 
PHI22 
PI180 


REAL *4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 


Description 


cos(Ls) coefficient 
constant coefficient of Al term 
latitude coefficient of Al term 
cos (2 Ls) coefficient 
constant coefficient of A2 term 
latitude coefficient of Al term 
square of latitude 
Ls phase angle of cos(Ls) term 
constant coefficient of PHI1 term 
latitude coefficient of PHI1 term 
Ls phase angle of cos (2 Ls) term 
constant coefficient of PHI2 term 
latitude coefficient of PHI2 term 
n / 180 


Mabhodology : 


Computes amplitudes of cos(Ls) and cos(2 Ls) terms as a function of 
the square of the latitude by 


al = all + al2*Lat2 
a2 = a21 + a22*Lat2 


PRSE 10 
PRSE 11 


Computes Ls phase angles of cos(Ls) and cos (2 Ls) terms by 


phil = phill + phil2 *Lat2 
phi2 = phi2 1 + phi22*Lat2 


PRSE 12 
PRSE 13 
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Compute the relative pressure, as a 


function of latitude and Ls - Lsun by 


PR =1. + al*cos (pil80* (Lsun - phil) ) + 
& a2*cos (2 . *pil80* (Lsun - phi2)) 


PRSE 14 
PRSE 15 
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Subroutine : Psur f ace 


Code: PSRF 


Description: Computes the daily average surface pressure and local surface 
pressure at a given latitude and longitude for a given areocentric 
longitude of the sun (Ls) , and daily average surface temperature. 

Called By: ATMOS 2 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

None 

Input Variables (not passed through commons): 


Name 

Type 

Description 

ALS . . 

. . REAL *4 

areocentric longitude of the sun (Ls) 

DUSTA . 

. . REAL* 4 

dust storm effect of diurnal amplitude 

DUSTM . 

. . REAL* 4 

dust storm effect on daily mean values 

GAMMA 

. . REAL *4 

local temperature lapse rate (K/km) 

GZERO . 

. . REAL *4 

surface gravity 

RBAR . . 

. . REAL* 4 

gas constant 

SI TLA . 

. . REAL *4 

local site latitude 

SITLO . 

. - REAL *4 

local site longitude 

SUNLON . 

. - REAL *4 

longitude of sub-solar point 

TAVG . . 

. . REAL* 4 

daily average surface temperature 

THGT . . 

. . REAL* 4 

local terrain height 

Output 

Variables (not 

passed through commons ) : 

Name 

Type 

Description 


PAVG. 

PSURF 


REAL *4 
REAL *4 


array of daily average pressures at significant levels 
daily average local surface pressure 


Local Variables (not passed through commons): 


Name 
All . . 

— 

Type 
REAL *4 

A12 . . 


REAL *4 

A21 . . 


REAL *4 

A22 . . 


REAL *4 

ABSLAT . 


REAL* 4 

ALSM60 . 


REAL *4 

AMP1 . . 


REAL *4 

AMP2 . . 


REAL *4 

CO thru 

C2 

REAL *4 

COSLAT . 


REAL *4 

COSLS . 


REAL *4 

DO thru 

D2 

REAL *4 

EXPONENT . . 

REAL *4 

FACTOR . 

. 

REAL *4 

FREQ . . 


REAL *4 

Hi 1 . . 


REAL *4 

H12 . . 


REAL *4 

H21 . . 


REAL *4 

H22 . . 


REAL *4 

NMALS . 


REAL *4 

NMLS. . 


REAL *4 

NPLS . . 


REAL *4 

PHI . . 

. 

REAL *4 


Description 


diurnal amplitude coefficient 
diurnal amplitude coefficient 
semi-diurnal amplitude coefficient 
semi-diurnal amplitude coefficient 
absolute value of latitude 
Ls - 60 degrees, converted to radians 
amplitude of diurnal dust effect 
amplitude of semi-diurnal dust effect 
coefficients in pressure/ temperature ratio 
cosine of latitude 
cosine of Ls - 60 

coefficients in pressure/ temperature ratio 

exponent in pressure variation with temperature 

coefficient in latitude and seasonal daily pressure 

diurnal frequency (15 degrees per Mars hour) 

diurnal phase coefficient 

diurnal phase coefficient 

semi-diurnal phase coefficient 

semi-diurnal phase coefficient 

square of 90 - absolute latitude 

square of 90 - latitude 

square of 90 + latitude 

time phase of diurnal surface pressure variation 
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PH2 . . 
PI180 . 
POTAVG . 
PR. . . 
PREF. . 
RREF . . 
SINLS . 
TBAR. . 
TLOCAL . 
TREF. . 
XLAT . . 


REALM 
REALM 
REALM 
REAL* 4 
REALM 
REAL *4 
REALM 
REALM 
REAL* 4 
REAL* 4 
REALM 


time phase of semi-diurnal surface pressure variation 
K / 180 

annual average pressure/ temperature ratio 

seasonally and latitudinally dependent relative pressure 

daily average temperature on reference ellipsoid 

local radius of reference ellipsoid 

sine of Ls 

estimated annual average temperature versus latitude 
local time in Mars hours (l/24th Sols) 

temperature in conversion of pressure to terrain height 
current latitude in radians 


Methodology: 


- 

latitude, using the RELLIPS subroutine. Evaluates the 
pressure / temperature ratio on the reference ellipsoid 


If (abslat .It. 55.) then 

potavg = 0 . 5* ( 2 . *C0 + Cl* (npls+nmls ) + C2*(npls**2 
& +nmls* *2 ) ) 

Else 

potavg = 0.5* (DO + Dl*abslat + D2*abslat**2 + 

Sc CO + Cl*nmals + C2*nmals**2) 

End if 


the 


PSRF 

29 

PSRF 

30 

PSRF 

31 

PSRF 

32 

PSRF 

33 

PSRF 

34 

PSRF 

35 


Computes the annual average temperature as a function of latitude 

ipbar = ( { ( 0 . 525161E-6*sitla - 0.263317E-5)*sitla 

& - 0 . 0116584) *sitla + 0 . 0334196 ) *sitla + 216.54 


PSRF 37 
PSRF 38 


Uses the PRSEAS subroutine to 
Ls value, then converts this 
reference ellipsoid 


get the relative pressure for the given latitude and 
to actual daily average pressure at the height of the 


Call PRSEAS (als , sitla, pr ) 


PSRF 40 


factor = 1.3854 - 3 . 59437E-5*sitla**2 
pref = f actor*pr*potavg*Tbar 


PSRF 43 
PSRF 44 


Converts pressure on the reference ellipsoid to pressure on the local terrain 
surface 


If (abs (gamma) .It. 0.001) then 

pavg = pref *exp (-1000. *gzero*thgt/ (Rbar*Tavg) ) 

Else 

Tref = Tavg + gamma* thgt 

exponent = -1000 . *gzero/ (Rbar* gamma) 

pavg = pref* ( (Tref /Tavg) **exponent) 

End if 


PSRF 46 
PSRF 47 
PSRF 48 
PSRF 49 
PSRF 50 
PSRF 51 
PSRF 52 


Computes the effects of dust storm (if any) on the diurnal average pressure 


pavg = pavg + ( 61 . 082 *sin( 2 . *pil80*sitla) 

& - 10 . 837*sin(4 . *pil80*sitla) ) *dustM 


PSRF 54 
PSRF 55 
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Evaluates the local time in Mars hours 


tlocal = 12. + (sunlon - sitlo)/15. 

If (tlocal .It. 0.) tlocal = tlocal + 24. 
If (tlocal .gt. 24.) tlocal = tlocal - 24. 


PSRF 57 
PSRF 58 
PSRF 59 


Computes the amplitude (ampl, N/m 2 ) 
in surface pressure 


and phase (phi, hours) of the diurnal variation 


sinLs = sin (pil80*als ) 

All = 12.60 - 0 . 118*sitla 
A12 = -0.60 - 0 . 063*sitla 
ampl = All + A12*sinLs 
Hll = 7.769 - 0 . 05769*sitla 
H12 = 3.077 + 0 . 001923 *sitla 
phi = Hll - H12 *sinLs 


PSRF 61 
PSRF 62 
PSRF 63 
PSRF 64 
PSRF 65 
PSRF 66 
PSRF 67 


Computes the amplitude (amp2, N/m**2) 
variation in surface pressure 


and phase (ph2, hours) of the semi-diurnal 


A21 = 7.982 - 0 . 022 *si tla 
A22 * -5.219 + 0 . 033 *si tla 
amp2 = A21 + A22*sinLs 
H21 = 11.635 - 0 . 02885*sitla 
H22 = -0.135 + 0 . 02885 *sitla 
alsm€0 = pil80* (als-60 . ) 
cosLs = cos(alsm60) 
ph2 = H21 + H22*cosLs 


PSRF 69 
PSRF 70 
PSRF 71 
PSRF 72 
PSRF 73 
PSRF 74 
PSRF 75 
PSRF 76 


Adds the dust storm effects on the diurnal and semi-diurnal amplitudes 


ampl = ampl + (1.32 + 13 . 56*coslat ) *dustA 

amp2 = amp 2 + (33 . 33*exp (- . 618727E-3 *sitla**2 ) ) *dustA 


PSRF 80 
PSRF 81 


Evaluates the local surface pressure as the daily average value 
semi-diurnal variations of pressure about the daily average 


plus the diurnal and 


psurf = pavg + ampl*cos (freq* (tlocal-phl) ) 
& + amp2*cos ( 2 . *freq* ( tlocal-ph2 ) ) 


PSRF 85 
PSRF 86 
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Funct i on : Random 


Cods: RAND 


Description s Returns a pseudo-random number, rectangularly distributed 
between 0 and 1 . 

Called By: MarsGRAM (Main) , Datastep, SETUP 

Connon Blocks Ussd (Sm Table 2-3 for a list of variables in coonon blocks): 

RANDCOM 

Input Variablss (not passsd through commons) : 

Name Type Description 

None 


Output Variablss (not passsd through conmons) : 


Name Type Description 

L INTEGER* 4 error flag 

RANDOM. . . REAL*4 uniformly distributed random variate (0 - 1) 


Local Variablss (not passsd through conmons) : 

Name Type Description 


ONE .... REAL* 4 1.0 

ZERO. . . . REAL *4 0.0 


Methodology: 

Returns a pseudo random variate, uniformly distributed between 0 and 
1, using the method of Algorithm AS 183, Applied Statistics (1982), Volume 
31, page 188. Uses the integer variables IX, IY and IZ, passed through 
common RANDCOM. Successive values of IX, IY and IZ are computed by 


IX = 171 * Mod(IX, 177) - 2 * (IX / 177) RAND 19 
IY = 172 * Mod(IY, 176) - 35 * (IY / 176) RAND 20 
IZ = 170 * Mod ( IZ , 178) - 63 * (IZ / 178) RAND 21 

If (IX .It. 0) IX = IX + 30269 RAND 23 
If (IY .It. 0) IY = IY + 30307 RAND 24 
If (IZ .It. 0) IZ = IZ + 30323 RAND 25 


The pseudo-random deviate value is computed from 


Random = Amod ( float ( IX) / 30269.0 + float(IY) / 30307.0 + RAND 37 

Sc float ( IZ) / 30323.0, one) RAN 13 38 
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Subroutine: RELLIPS „ _ 

Code : RLPS 

Description: Calculates the reference radius of Cain et al's 6 1 mbar 
reference ellipsoid, for a given latitude. Also computes the acceleration 
latitude ^ or a 9 lven altitude above the reference ellipsoid at that 

SETO* d By! MarSGRAM (Main) ' ATM0S2 ' Datastep, Psurface, ST EWART 2 , THERMOS, 


Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

None 

Input Variables (not passed through commons) : 

^ ame Type Description 


LAT . . 

Z . . . 

. . REAL *4 

• . REAL* 4 

current latitude (degrees) 
current height (km) 


Output 

Variablos (not passed through commons) * 


Name 

Type 

Description 


GZ. . . 

RREF . . 

. . REAL *4 

. . REAL *4 

acceleration of gravity at height z 
local radius of reference ellipsoid 

(m/s 2 ) 

(km) 

Local Variables (not 

passed through commons) : 


Name 

Type 

Description 



A . 
AB. 

B . 

C . 
GM. 

J2 . 

P2 . 
RADEG 
RZ . 
TLAT 
XX. 
YY. 


REAL *4 
REAL *4 
REAL *4 
REAL * 4 
REAL *4 
REAL* 4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 


largest equatorial radius of reference ellipsoid 
square root of product of A and B 
smallest equatorial radius of reference ellipsoid 
polar radius of reference ellipsoid 
gravitational constant time mass of Mars 
coefficient of 1st non-spherical component of gravity 
latitude part of non-spherical component of gravity 
180 / n 

radius to current height 
tangent of latitude 

square of x component of local ellipsoid radius 
square of y component of local ellipsoid radius 


Methodology : 

3393 21° and in r t - y ^ 7 fi a 7 B e l liPSOid * nd the reference ra <*:Li A = 3394.67 km, B = 
3393.21 km, and C - 3376.78 to compute the local radius of the reference ellipsoid 


XX - (AB * C)**2 / (C**2 + (AB * TLAT)* *2) 
YY = XX * TLAT* *2 
Rref = SQRT (XX + YY) 


RLPS 15 
RLPS 16 
RLPS 17 


Compute the acceleration of gravity at altitude 
in Space Research, Volume 2, 1982, pages 3-17). 


z, by the method of Seiff (Advances 
This includes the spherically 
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symmetric gravitational component and the first <J2) non-spherically symmetric 
component 


GM = 4 . 28282E7 
J2 = 0.001965 

P2 = 1 . 5 * s in ( LAT / RADEG ) * * 2 - 0.5 
Rz = Rref + z 

gz = (GM/Rz**2)Ml- - 3 . * J2 M (AB/Rz) * *2 ) *P2 ) 


RLPS 19 
RLPS 20 
RLPS 21 
RLPS 22 
RLPS 23 
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Subroutine s St ratos 


Codex STRA 


Description: Interpolates between z75, the areographic height for the 75 km 
areopotential significant level, and ZF, the base of the Stewart 
thermosphere. Computes temperature, pressure, density, scale height, and 
temperature gradient at a given height z between these two levels. 

Callsd By: ATMOS 2 

Common Blocks Used (See Tabls 2-3 for a list of variablss in common blocks) : 

None 

Input Variablss (not passsd through commons) : 

Name Type Description 


acceleration of gravity at current altitude 
molecular weight at 75 km areopotential km 
molecular weight at base of thermosphere 
pressure at 75 km areopotential km 
pressure at base of thermosphere 
local radius of reference ellipsoid 
gas constant 

temperature at 50 km areopotential km 
temperature at 75 km areopotential km 
temperature at base of thermosphere 
current altitude (km) 

areographic height at 50 km areopotential km 
areographic height at 75 km areopotential km 
height at base of thermosphere 


Output Variablss (not passsd through commons) : 

Name Type Description 


. . REALM atmospheric density at current altitude 

. . REAL* 4 scale height at current altitude 

. . REAL* 4 atmospheric pressure at current altitude 

. . REAL *4 vertical temperature gradient (K/km) 

* . REAL* 4 atmospheric temperature at current altitude 


Local Variablss (not passsd through commons) : 

Name Type Description 


height interpolation factor (H - H75)/(HF - H75) 

exponent in pressure versus temperature dependence 

temperature gradient between H75 and HM 

temperature gradient between HM and HF 

gravity divided by gas constant 

areopotential altitude corresponding to Z 

areopotential altitude corresponding to Z50 

areopotential altitude corresponding to Z75 

areopotential altitude corresponding to ZF 

HF adjusted by molecular weight ratio 

midpoint between H75 and HFP 

H adjusted by molecular weight ratio 

molecular weight at height Z 

pressure at midpoint height HM 

temperature at midpoint height HM 



HZ. . 
PZ. . 
TGRAD 
TZ. . 


M75 . 
MF. . 
P75 . 
PF. . 
RREF . 
RSTAR 
T50 . 
T75 . 
TF. . 
Z . . 
Z50 . 
Z75 . 
ZF. . 


REAL* 4 
REAL* 4 
REAL* 4 
REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL* 4 
REAL* 4 
REAL* 4 
REAL*4 
REALM 
REAL *4 
REAL *4 
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Methodology i 


Uses input values of 275 and z50, the areographic altitudes above the 
reference ellipsoid corresponding to 75 km and 50 km areopotential altitude above 
local surface, to compute H75 and H50, the areopotential altitudes above the 
reference ellipsoid, corresponding to z75 and Z50 


H75 = z75/(l. + z75/rref ) 
H50 = z50/(l. + Z50/rref ) 


STRA 8 
STRA 9 


Finds HF, the areopotential altitude corresponding to zF, areographic altitude of 
the base of the thermosphere 


HF = zF/(l. + zF/rref) 


STRA 12 


Computes H, the areopotential altitude corresponding to the current areographic 
altitude z, relative to the reference ellipsoid 


H = z/ (1 . + z/rref ) 


STRA 15 


Interpolates on height to get M, the molecular weight at height z 


dH = (H - H75 ) / (HF - H75) 
M = M75 + (MF - M75 ) *dH 


STRA 17 
STRA 18 


Adjusts the heights HF and H by the molecular weight ratio 


HFP = MF*HF/M75 
HP = M*H/M75 


STRA 20 
STRA 21 


Finds HM, the midpoint (molecular weight adjusted) height between z75 and ZF 


HM = 0.5* (HFP + H75) 


STRA 23 


Uses quadratic height interpolation (fitting the temperature at heights H50, H75 and 
HFP) to get the temperature TM at the midpoint height HM 


TM = -0 . 25*T50* (HFP-H75) **2/ ( (H75-H50) * (HFP-H50) ) 

& + 0 . 5*T75* (HM-H50 ) / (H75-H50 ) + 0 . 5*TF* (HM-H50) / (HFP-H50) 

If (TM. eq . T75 ) TM = T75 + 1. 

If (TM. eq . TF) TM = TF - 1 . 


STRA 25 
STRA 26 
STRA 27 
STRA 28 


Computes the temperature gradients for the two layers H75 to HM and HM to HF 


garni = (TM - T75)/(HM - H75) 
garn2 = (TF - TM) / (HFP - HM) 


STRA 30 
STRA 31 


Find the exponent factors for pressure variation in the two layers 
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goR = Alog (p75/pF) / (Alog (TM/T75) /gaml + Alog (TF/TM) /gam2 ) 
pM = p7 5 * (T75/TM) * * (goR/gaml ) 


STRA 33 
STRA 34 


Interpolates for temperature Tz and pressure pz at height z, using linear 
temperature gradients in whichever of the two layers the height z falls (H75 to HM 
or HM to HFP ) 


If (HP .It. HM) Then 
tgrad = gaml 

Tz = T75 + tgrad* (HP - H75) 
expon = goR/gaml 
pz = p7 5 * (T75/Tz ) * *expon 
Else 

tgrad = gam2 

Tz = TM + tgrad* (HP - HM) 
expon = goR/gam2 
pz = pM* (TM/Tz) **expon 
Endif 


STRA 37 
STRA 3 8 
STRA 39 
STRA 40 
STRA 41 
STRA 42 
STRA 43 
STRA 44 
STRA 45 
STRA 46 
STRA 47 


Convert the units of the temperature gradient to K/km 


tgrad = 0.001*tgrad 

Computes the density dz at height z from the perfect gas law 
100 dz = M*pz/ (Rstar*Tz) 

Calculates the scale height Hz at height z by the relation 
Hz = 0 . 001 *Rs tar *Tz/ (gz*M) 


STRA 49 


STRA 51 


STRA 53 
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Subroutine : STEWART2 


Code: STW2 


Description: Stewart thermosphere model for computing global mean 
temperature and density at a given areocentric longitude of the sun (Ls) . 
Computes effects due to solar activity influence and effects due to dust 
storm (if one is in progress) . 

Called By: ATMOS 2 

Common Blocks Used (See Table 2-3 Cor a list of variables in common blocks): 


THERM 


Input Variables (not passed through commons): 


Name Type 


Description 


ALSO. 
CHGT . 
INTENS 
IU0 . 
LAT . 
LST . 
LSUN . 
RAUI . 
RSTAR 
SIGMA 


REAL *4 
REAL* 4 
REAL* 4 
INTEGER* 4 
REAL* 4 
REAL* 4 
REAL *4 
REAL* 4 
REAL* 4 
REAL *4 


Ls value at start of dust storm (if any) 
current height 

dust storm intensity (0.0 - 3.0) 
unit number for screen output of messages 
current latitude 
local solar time 

areocentric longitude of the sun (Ls) 

Mars orbital radius in astronomical units 
gas constant 

standard deviation for long term thermosphere variations 


Output Variables (not passed through commons) : 


Name 

Type 

Description 

H . . . . 

REAL *4 

pressure scale height 

MOLWTG . . 

REAL* 4 

molecular weight 

TOTALMDZ. 

REAL* 4 

total atmospheric mass density 

TOTALPRZ . 

REAL* 4 

total atmospheric pressure 

TZ .... 

REAL*4 

atmospheric temperature 


Local Variables (not passed through commons) : 


Name Type 


Description 


DR. . 

. 

. * 

REAL* 4 

DUST. 

• 

* • 

REAL*4 

ES. . 



REAL *4 

FBAR . 



REAL* 4 

FBARR 



REAL *4 

FLAG. 



INTEGER* 4 

GZ. . 



REAL* 4 

MDZ . 



REAL* 4 

NDZ . 



REAL* 4 

PFAC. 



REAL *4 

PRZ . 



REAL* 4 

RAU . 



REAL* 4 

RF. . 



REAL* 4 

RREF . 



REAL* 4 

SMA . 



REAL* 4 

TF. . 



REAL* 4 

TINF. 



REAL *4 

TO. . 

. 


REAL* 4 


seasonal correction to the height of the base of the 
thermosphere 

dust correction to the height of the base of the 
thermosphere 

array of thermospheric correction factors 

mean value of 10.7 cm solar flux at Earth position (1 AU) 
mean value of 10.7 cm solar flux at Mars position 
test flag, set to 1 to get diagnostic output 
acceleration of gravity at current height 
atmospheric mass density 
atmospheric number density 

relative pressure factor for season and latitude 

atmospheric pressure 

Mars orbital radius (in AU) 

total radius to position (height plus ellipsoid radius) 

local radius of reference ellipsoid 

semi-major axis of Mars orbit 

temperature at the base of the thermosphere 

exospheric temperature 

factor used in correction to ZF for seasonal variation 
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REAL *4 
REAL *4 
REAL* 4 


total atmospheric number density 
height of the base of the thermosphere 
difference between current height and ZF 


TOTALNDZ . . 

ZF 

ZZF . . . . 


Methodology : 

The Stewart thermosphere model (including subroutine STEWART 2 and the 
associated subroutines THERMOS , DZDUST, Escalc, PRSEAS and RELLIPS ) was converted to 
FORTRAN code from the Pascal version listed in Appendix B of "The Mars Atmosphere: 
Observations and Model Profiles for Mars Missions*, David E. Pitts et al., eds 
JSC-24455. The model has also been discussed, and given in IDL code, by Ian 
Stewart, Laboratory for Atmospheric and Space Physics, University of Colorado Final 
Report JPL PO # NQ-802429. 

Use the EScalc subroutine to evaluate the ES array of thermospheric 
variability factors for the input standard deviations for short term (stdl) and long 
term (SIGMA) variability 


Call EScalc (stdl , SIGMA, ES) 


STW2 21 


Adjust the 10.7 cm solar flux for long-term mean value, with ES(0) 


FBAR = F107 * EXP (ES ( 0 ) ) 


STW2 25 


Convert the 10.7 cm solar flux to the orbital position of Mars 


FBARR = FBAR / (RAU**2) 


STW2 29 


Call RELLIPS to get the acceleration of gravity, GZ, and the local radius of the 
reference ellipsoid, RREF . Call PRSEAS to get the relative pressure factor for the 
given value of Ls (LSUN) and latitude. 


CALL RELLIPS ( LAT, RREF, CHGT, GZ) STW2 30 

CALL PRSEAS (LSUN, LAT, PFAC) STW 2 31 


Calculate DR, the seasonal correction to the height of the base of the thermosphere 


TO = 220.0 * SMA / RAU 

DR = (TO / 19.51) * ALOG(PFAC) 


STW2 33 
STW2 34 


Evaluate the dust storm correction and compute ZF, the height of the base of the 
thermosphere, and ZZF, the current height relative to the base of the thermosphere 


CALL DZDUST ( LSUN, alsO, INTENS, DUST) 

DUST = DUST * EXP (ES (10) ) 

ZF = (124.4 * SMA / RAU) * EXP(ES{8) + ES ( 9 ) ) + DR + DUST 
ZZF = CHGT - ZF 


STW2 36 
STW2 37 
STW2 38 
STW2 39 


Find the total radius (radius of reference ellipsoid plus the height), the 
exospheric temperature, TINF, and the temperature of the base of the thermosphere 
TF 


RF - RREF + ZF STW2 45 

TINF = 4.11 * (11.0 + FBARR) * EXP(ES(2) + ES(3)) STW2 46 
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TF = (170.0 * SMA / RAU) * EXP(ES(8) + ES(9)) 


STW2 47 


Call the THERMOS subroutine to get the remaining variables computed by 
model 


CALL THERMOS (FLAG, ES, TINF, TF, LAT, LST , ZF , RF , ZZF, TOTALPRZ , 
St TOTALNDZ, TZ, MOLWTG , PRZ, NDZ , MDZ , TOTALMDZ ,iu0) 


Calculate H, the scale height, in km 


H = RSTAR*TZ/ ( 1000 . *MOLWTG*GZ) 


Convert the atmospheric pressure value to N/m**2 


TOTALPRZ = TOTALPRZ*! . 0E5 


Convert the atmospheric density value to kg/m**3 


TOTALMDZ = TOTALMDZ * 1 0 0 0 . 


the Stewart 


STW2 58 
STW2 59 


STW2 61 


STW2 63 


STW2 65 
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Function: 


Tdiurnal 


Code: TDIR 


Description: Computes local surface temperature from the daily maximum and 
minimum surface temperature and the site-sun longitude separation. 

Called By: Tsurface 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

None 


Input Variables (not passed through commons) : 

Description 


Name 

Type 

DLO . . 

. . REAL* 4 

TMAX. . 

. . REAL *4 

TMIN. . 

. . REAL* 4 

Output 

Variables (not 

Name 

Type 

TDIURNAL. , REAL * 4 

Local Variables (not 

Name 

Type 

COSFAC . 

. . REAL *4 

DT. . . 

. . REAL *4 

OFFSET . 

. . REAL *4 

PERIOD. 

. . REAL *4 

PI180 . 

. . REAL *4 


longitude difference from current position to sun 
daily maximum surface temperature 
daily minimum surface temperature 


Description 

local surface temperature at current time 


Description 


cosine factor in diurnal variation (+1 or -1) 
diurnal temperature interpolation factor 
longitude offset in diurnal temperature interpolation 
effective period for diurnal interpolation 
7t / 180 


Methodology : 


Divides site-sun longitude difference into five ranges (-180 to -102.19 
degrees, -102.19 to -48.6 degrees, -48.6 to +43.8 degrees, +43.8 to +111.5 degrees, 
and +111.5 to +180 degrees). Finds values of cosfac, period and offset for 
whichever range is appropriate, based on input value of site-sun longitude 
difference. Computes the diurnal temperature shape factor dT for the values of 
cosfac, period and offset thus determined 


dT = . 5*(1. + cosfac*cos <pil80*period* (dlo + offset))) TDIR 26 

dT = {-0.3455*dT + 1.3455) *dT TDIR 27 


Evaluates the local surface temperature, from the interpolation factor dT and the 
daily maximum and minimum temperature values 


Tdiurnal = Tmin + (Tmax - Tmin) *dT 


TDIR 29 
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Functions 


Terrain 


Code : TERN 


Description: Computes terrain height, relative to reference ellipsoid, for 
the given array of terrain height data, at a given latitude and longitude. 

Called By: MarsGRAM (Main), Datastep, SETUP 

Common Blocks Used (See Table 2-3 for a list of variables in common blocks) : 

TERHGT 

Input Variables (not passed through commons) : 


Name Type Description 


ALAT. . . 

REAL* 4 

current 

latitude ( degrees ) 

ALON . . . 

REAL*4 

current 

West longitude (degrees) 


Output 

Variables (not 

passsd through commons): 

Name 

Type 

Description 

TERRAIN 

. . REAL*4 

terrain height (km) 


Local Variables (not paaaad through common* ) : 

Name Type Description 

DLAT . . 


REAL*4 

latitude increment between terrain data base values 

DLON . . 


REAL*4 

longitude increment between terrain data base values 

FI thru 

F4 

REAL*4 

terrain heights at four corners of grid containing 
the desired location 

ILAT . . 


INTEGER* 4 

latitude index of terrain height data array 

ILON . . 


INTEGER* 4 

longitude index of terrain height data array 

NLAT. . 


INTEGER* 4 

number of latitude grids in data array 

NLON. . 


INTEGER* 4 

number of longitude grids in data array 

RADVL1 . 


REAL* 4 

normalized distance parameter from Viking 1 site 

RADVL2 . 


REAL* 4 

normalized distance parameter from Viking 2 site 

TERR. . 


REAL *4 

intermediate value of interpolated terrain height 

VL1HGT . 


REAL* 4 

Viking 1 site terrain height 

VL1LAT. 


REAL* 4 

Viking 1 site latitude 

VL1LON . 


REAL* 4 

Viking 1 site longitude 

VL2HGT . 


REAL *4 

Viking 2 site terrain height 

VL2LAT . 


REAL*4 

Viking 2 site latitude 

VL2LON. 


REAL *4 

Viking 2 site longitude 

XLON . . 


REAL *4 

relative longitude from interpolation grid corner 

YLAT. . 


REAL* 4 

relative latitude from interpolation grid corner 


Methodology : 

Computes the latitude and longitude increments between the terrain data bases 
grid points 


dlat = 180 . / {nlat-1 * ) TERN 14 

dlon = 3 60 . / (nlon-1 . ) TERN 15 


Finds the latitude and longitude index values of the interpolation lat-lon rectangle 
that contains the desired location 
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ilat = 1 + int ( (90.+alat) /dlat) 

If (ilat .gt. nlat-1 ) ilat = nlat-1 


TERN 17 
TERN 18 


ilon = 1 + int (alon/dlon) 

If (ilon .gt. nlon-1) ilon = nlon-1 


TERN 20 
TERN 21 


Gets the terrain heights at the corner points of the interpolation rectangle from 
the terrain height data base array 


FI = th { ilat , ilon) 

F2 = th ( ilat , ilon+1 ) 

F3 = th ( ilat+1 , ilon) 

F4 = th ( ilat+1 , ilon+1 ) 


TERN 23 
TERN 24 
TERN 25 
TERN 26 


Computes normalized latitude and longitude displacements from the corner of the 
interpolation rectangle 


ylat = 1. + (90.+alat)/dlat - ilat 
xlon = 1. + (alon/dlon) - ilon 


TERN 28 
TERN 30 


Does a bi-linear interpolation across the interpolation rectangle to find the 
terrain height at the desired location 


Terr = FI + (F2 - FI) ♦xlon + <F3 - FI) *ylat 
& + (F4 - F2 - F3 + Fl ) *xlon*ylat 


TERN 32 
TERN 33 


Finds the normalized distances from the Viking 1 and Viking 2 sites 


radvll — ( ( alat- vl Hat ) * * 2 + (alon-vlllon) **2 ) /4 . 0 TERN 35 

radvl2 = { (alat-vl21at ) **2 + (alon-vl21on) **2 ) /4 . 0 TERN 36 


Modifies the interpolated terrain value so that there is a smooth transition to the 
Viking site value, if the location is within about 1 degree of lat-lon from either 
of the Viking sites 


If (radvll .It. 1.0) then 

Terr = { 1 . -radvll ) *vllhgt + Terr*radvll 
Else If(radvl2 .It. 1.0)then 

Terr = ( 1 . -radvl2 ) *vl2hgt + Terr*radvl2 
Endif 

Terrain = Terr 


TERN 38 
TERN 39 
TERN 40 
TERN 41 
TERN 42 
TERN 43 
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Subroutine t THERMOS Code: THRM 

Description: Subroutine for the Stewart thermosphere model. Returns the 
temperature, pressure, molecular weight, mass density and number density 
concentrations versus altitude above the base of the thermosphere . 

Callsd By: ST EWART 2 

Common Blocks us®d (Sss Table 2-3 for a list of variables in c numn n blocks) : 

None 


Input Variables (not passed through commons) : 


Name 


Type 

Description 



ES .... 


REAL* 4 

array of thermospheric variability factors 

for 

the input 




standard deviations for short term (stdl) 

and 

long term 




(SIGMA) variability 



FLAG . . . 


REAL* 4 

test flag, set to 1 to get diagnostic output 


IUO . . . 


INTEGER* 4 

unit number for screen output of messages 



LAT . . . 


REAL* 4 

current latitude 



LST . . . 


REAL *4 

local solar time 



RF. . . . 


REAL*4 

total radius (radius of reference ellipsoid plus height) 

TF .... 


REAL*4 

temperature at base of thermosphere 



TINF . . . 


REAL *4 

exospheric temperature 



ZF .... 


REAL *4 

height of base of thermosphere 



ZZF . . . 

• 

REAL *4 

current height above height ZF 



Output Variables (not 

passed through cannons): 



Name 


Type 

Description 



MDZ . . . 


REAL* 4 

atmospheric constituent mass density 



MOLWTG . . 


REAL *4 

molecular weight 



NDZ . . . 


REAL *4 

atmospheric constituent number density 



PRZ . . . 


REAL *4 

atmospheric constituent partial pressure 



TOTALMDZ. 


REAL *4 

total atmospheric mass density 



TOTALNDZ . 


REAL* 4 

total atmospheric number density 



TOTAL PRZ . 


REAL* 4 

total atmospheric pressure 



TZ .... 


REAL *4 

local temperature at current height 




Local 

Name 

Variables (not passad through coanons) : 

Type Description 

AO. . 



REAL *4 

factor in atomic oxygen calculation 

BK. . 



REAL *4 

numerical constant for the number density calculations 

DM. . 



REAL *4 

molecular mass array 

FF. . 



REAL *4 

fractional composition of the heterosphere 

FO. . 



REAL* 4 

factor in atomic oxygen calculation 

GF. . 



REAL* 4 

acceleration of gravity at height ZF 

HH. . 



REAL *4 

scale heights of the heavy gas constituents 

I . . 



INTEGER* 4 

index for heavy gas calculations 

J . . 



INTEGER*4 

index for light gas calculations 

K . . 



INTEGER* 4 

dummy index in implied loop for write statement 

M . . 



REAL *4 

molecular weights for the constituent gases 

Pi BAR 



REAL *4 

pressure of 1 bar 

PFH . 



REAL* 4 

partial pressure of H at height ZF 

PFH2 . 



REAL* 4 

partial pressure of H2 at height ZF 

PFHE . 



REAL*4 

partial pressure of He at height ZF 

PRESSF . 

. . 

REAL* 4 

pressure at height ZF 
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RADEG 
RATIO 
RREF . 
SCALE 
XDM . 
XFF . 
XHH . 
YSC . 


REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 


180 / n 

factor used in computing partial pressure of H 

local radius of reference ellipsoid 

scale factor used in temperature calculation 

molecular mass array for light gases 

partial pressure array for He, H2 and H 

scale heights of the light gas constituents 

scaled height parameter in the temperature calculation 


Methodology : 


The Stewart thermosphere model (including the main subroutine STEWART 2 and 
the associated subroutines THERMOS , DZDUST, Escalc, PRSEAS and RELLIPS) was 
converted to FORTRAN code from the Pascal version listed in Appendix B of "The Mars 
Atmosphere: Observations and Model Profiles for Mars Missions", David E. Pitts et 
a * ' eds . , JSC-24455. The model has also been discussed, and given in IDL code by 
Ian Stewart, Laboratory for Atmospheric and Space Physics, University of Colorado, 
Final Report JPL PO # NQ-802429. 

Calls RELLIPS to get local radius of reference ellipsoid, RREF, and GF the 
acceleration of gravity at the height ZF, the base of the thermosphere 


CALL RELLIPS (LAT, RREF, ZF,GF) 


THRM 27 


Sets the pressure at the base of the thermosphere to 1.26 nanobars 
PRESSF = 1.26E-9 

Sets PI BAR to the pressure of 1 bar (1.0x10 s N/m 2 ) 

P1BAR = 1 . 0E6 


(1.26X1CT 4 N/m 2 ) 
THRM 31 


THRM 34 


Computes AO and F0, parameters for use in calculating the atomic oxygen 
concentration 


AO = 0.18 * (1.0 + ES ( 7 ) ) 

FO = 0.01 * EXP ( ES ( 4 ) + ES ( 5 ) ) 


THRM 3 6 
THRM 3 8 


Evaluates SCALE, the scale factor used in the temperature versus height 


SCALE = TF / 9 . 2 0 


THRM 40 


Sets the molecular weight constant arrays M, DM and XDM (THRM 41 through THRM 76). 
Sets the array FF, for the composition of the heterosphere (THRM 77 through THRM 
89) . Adjusts FF ( 5 ) for atomic oxygen 


FF(5) = FO* ( 1 . 0- AO* SIN ( 15 . 0*LST/ RADEG) *COS (LAT/ RADEG) ) 


THRM 90 


® xo ^ ase ( ^ e ^ ght ZF) Partial pressures for helium (PFHE) , molecular hydrogen 
PFH2 and atomic hydrogen ( PFH) and converts these to partial pressures (XFF array) 


PFHE = 3.3E-16 * TINF 


THRM 91 



THRM 93 


PFH2 = 2.4E-15 

if (TINF .LE. 330.0) PFH=5 . 2E-16*TINF*EXP (-TINF/70 . 0) 

if (TINF .GT. 330.0) then 
RATIO = 1440.0 / TINF 

PFH=5.BE-18*SQRT(TINF) *EXP(RATIO) / (l.O+RATIO) 


ENDIF 

XFF(O) = PFHE / PRESSF 
XFF(l) = PFH2 / PRESSF 
XFF ( 2 ) = PFH / PRESSF 


THRM 95 

THRM 97 
THRM 98 
THRM 99 

THRM101 

THRM102 

THRM103 

THRM104 


Calculates local atmospheric temperature at current height 


YSC = ZZF * RF / (RF + ZZF) 

TZ = TINF - (TINF - TF) * EXP(-YSC / SCALE) 


THRM108 

THRM109 


Loops through the heavy gases (I = 0 for C02 , 1 for N2. 2 for Argon 3 for 02, 4 for 
CO and 5 for atomic oxygen), and computes scale height, HH(I), partial pressure, 
PRZ(I), number density, NDZ(I), mass density, MDZ(I), and accumulates total mass 
density, TOTALMDZ , total pressure, TOTALPRZ , and pressure-weighted molecular weight, 

MOLWTG. 


DO 200 I = 0, 5 

HH ( I ) = BK * TINF / (GF * DM ( I ) ) / 1 . 0ES 


THRM110 

THRM111 


PRZ(I) = PRESSF’FFd) *EXP ( -YSC/HH < I ) - (SCALE/HH(I) ) *ALOG(TZ/TF) ) THRM113 
NDZ(I) = P1BAR * PRZ(I) / (BK * TZ) THRM114 


MDZ(I) = NDZ(I) * DM (I) 

TOTALMDZ = TOTALMDZ + MDZ ( I ) 
TOTALPRZ = TOTALPRZ + PRZ ( I ) 
MOLWTG = MOLWTG + PRZ ( I ) * M ( I ) 
200 CONTINUE 


THRM116 

THRM117 

THRM118 

THRM119 

THRM120 


Loops through the light gases (J = 0 for He, 1 for H2, 2 for H) and. scale 
heiaht XHH(J) , partial pressure, PRZ(J+6), number density, NDZ(J+6), mass density, 
MDZ^J+6 ) ^and accumulates total mass density, TOTALMDZ, total pressure, TOTALPRZ, 
and pressure-weighted molecular weight, MOLWTG. 


DO 210 J = 0, 2 

XHH(J) = BK * TINF / (GF * XDM(J)) / 1 . 0E5 


THRM121 

THRM122 


PRZ (J+6 ) =PRESSF*XFF ( J) *EXP ( -YSC/XHH ( J) - ( SCALE /XHH ( J) ) * 
& ALOG (TZ/TF) ) 

NDZ ( J+6 ) = P1BAR * PRZ ( J +6) / (BK * TZ) 


THRM124 

THRM125 

THRM126 


MDZ ( J +6) = NDZ ( J +6) * XDM(J) 
TOTALMDZ = TOTALMDZ + MDZ(J + 6) 
TOTALPRZ = TOTALPRZ + PRZ(J + 6) 

MOLWTG = MOLWTG + PRZ ( J +6) * M(J + 6) 
210 CONTINUE 


THRM128 

THRM129 

THRM130 

THRM131 

THRM132 


Divides by total pressure to get true total molecular weight, and computes total 
number density, TOTALNDZ 


MOLWTG = MOLWTG / TOTALPRZ 

TOTALNDZ = P1BAR * TOTALPRZ / (BK * TZ) 


THRM133 

THRM134 
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Subroutine: Temps „ . 

Code : TMPS 

©ascription: Computes array of atmospheric temperatures at the significant 
levels for a given latitude and longitude. Dust effects are included. 

Celled By: ATMOS 2 

Common Blocks Used (See Table 2-3 for a list of variables in o 

None 


blocks) : 


Input Variables (not passed through 


)i 


Name 


Type 


Description 


ALAT. . 

. . REAL* 4 

ALS . . 

. . REAL* 4 

DLON . . 

. . REAL *4 

DUSTA . 

. - REAL *4 

DU STM . 

. . REAL* 4 

GAM . . 

. . REAL* 4 

T0BAR . 

. . REAL *4 

TSURF . 

. . REAL* 4 

Output 

Variables (no 

Name 

Type 

T . . . 

- . REAL *4 

GAM . . 

. . REAL *4 

Local Variables (not 

Name 

Type 

A25 . . 

. . REAL *4 

ABSLAT. 

. . REAL *4 

AFACT . 

. . REAL* 4 

ALATP . 

. . REAL* 4 

ALONMAX 

. . REAL *4 

AT25. . 

. - REAL* 4 

B25 . . 

. . REAL *4 

BT2 5 . . 

. . REAL* 4 

C25 . . 

. . REAL *4 

COSLAT . 

. . REAL *4 

CT25 . . 

. . REAL *4 

DFACTOR 

. . REAL *4 

DT25 . . 

. . REAL *4 

DTR26 . 

. . REAL *4 

DTR30 . 

. . REAL *4 

DZ . 

- . REAL *4 

FACTOR . 

. . REAL *4 

I . . . 

. . INTEGER* 4 

PI120 . 

. . REAL *4 

PI180 . 

. . REAL *4 

SINLAT. 

. . REAL *4 

T25 . . 

. . REAL *4 

THOUR . 

. . REAL *4 

THOURO . 

. . REAL *4 

THOUR26 

. . REAL *4 

Z . . . 

. . REAL *4 


local site latitude 

areocentric longitude of the sun (Ls) 

site-sun West longitude difference 

magnitude of dust storm effect on daily amplitude 

magnitude of dust storm effect on daily average 

lapse rates (-dT/dz) between significant levels 

average surface temperature with no dust storm effect 

local surface temperature, including dust effect 


Description 


array of air temperatures at the significant levels 
lapse rates, modified by dust storm effects 


saoa tnrough commons) : 

Description 


coefficient array in calculation of T25 
absolute value of latitude 

factor array for diurnal temperature variation 

90 minus absolute latitude 

longitude where maximum phase occurs 

latitude dependent coefficient in T25 calculation 

coefficient array in calculation of T25 

latitude dependent coefficient in T25 calculation 

coefficient array in calculation of T25 

cosine latitude term in T25 calculation 

latitude dependent coefficient in T25 calculation 

factor array for dust storm effect on diurnal amplitude 

difference between T25 and value from lapse rates 

diurnal amplitude factor for 26 km 

diurnal amplitude factor for 30 km 

array of height differences between significant levels 
height scaling factors for diurnal amplitude (relative 
to amplitude at 30 km) in non-dust storm conditions 
array index 

7i / 120 
n / 180 

sine latitude term in T25 calculation 
air temperature at 25 km 

phase hour of maximum diurnal perturbation 

phase hour for maximum diurnal perturbation at surface 

phase hour for maximum at 26 km height 

height array for significant levels 


A-63 



Methodology: 


Starts by computing an estimate of the temperature at 25 km, T25, from a set 
of coefficients for latitude and Ls dependence 


AT25 = A25 (0) 

BT25 = B25 ( 0 ) 

CT25 = C2 5(0) 

Do 5 i = 1,5,2 

sinlat = sin (pil20* ( i+1 . ) *alat/2 . ) 
coslat = cos (pil20* ( i+1 . ) *alat/2 . ) 

AT25 = AT25 + A25 ( i) *sinlat 
AT25 = AT25 + A25 ( i+1) *coslat 
BT25 = BT25 + B25 ( i ) *sinlat 
BT25 = BT25 + B25 ( i+1 ) *coslat 
CT25 = CT25 + C25 (i) *sinlat 
5 CT25 = CT2 5 + C25 (i+1) *coslat 

T25 = AT25 + BT25*sin (pil80*als) + CT25*cos (pil80*als) 

Computes the daily average atmospheric temperature (no dust storm correction) at 
each of the significant levels, from the daily average surface temperature and the 
lapse rates between the significant levels 


TMPS 41 
TMPS 42 
TMPS 43 
TMPS 44 
TMPS 45 
TMPS 46 
TMPS 47 
TMPS 48 
TMPS 49 
TMPS 50 
TMPS 51 
TMPS 52 
TMPS 53 


T(0) = TObar 
Do 10 i = 1,3 

10 T(i) = T(i-l) - gam(i) *dz (i) 


TMPS 56 
TMPS 57 
TMPS 58 


Finds the difference from the temperature estimate T25 and the temperature at 25 km, 
computed from the lapse rates 

DT25 = T25 - (T ( 2 ) + 2.*T(3))/3. ™ ps 60 


Adjusts the temperature and lapse rates for the first three significant levels, 
based on the value of T25 


Do 15 i = 1,3 

T(i) = T(i) + z (i) *DT25/25 . 
15 gam(i) = (T(i-l) - T(i))/dz(i) 


value of T ( 3 ) 


Do 20 i = 4,5 

20 T ( i ) = T(i-l) - gam(i) *dz (i) 


km height 

If (dustM .eq. 0)then 

If (abslat .le. 25) then 
Factor(l) = 0.79 
Else If (abslat .ge. 45) then 
Factor ( 1 ) = 0.43 




TMPS 

62 



TMPS 

63 



TMPS 

64 

5th 

significant 

levels, based on the 

new 



TMPS 

66 



TMPS 

67 

= 0 

(TMPS 70 - 

TMPS 81) . Calculates 

factor 

for 

the diurnal 

temperature variation at 5 



TMPS 

70 



TMPS 

72 



TMPS 

73 



TMPS 

74 



TMPS 

75 
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Else 

Factor { 1 ) = 1.24 - 0.018*abslat 
Endif 

Do 25 i = 1,5 
25 afact(i) = factor(i) 

Endif 


TMPS 7 6 
TMPS 77 
TMPS 78 
TMPS 79 
TMPS 80 
TMPS 81 


Computes the relative temperature amplitude at 30 km height 


dtr3 0 = 0.0141* (90 . - abslat)**0.2 


TMPS 84 


Evaluates thour26, the phase hour for maximum temperature at 26 km height 


If (alat .le. -70.) then 
thour26 = 21. 

Else If (alat .ge. -20.)then 
thour26 = 11.5 
Else 

thour2 6 = 7.7 - 0.19*alat 
Endif 


TMPS 89 
TMPS 90 
TMPS 91 
TMPS 92 
TMPS 93 
TMPS 94 
TMPS 95 


Computes the effects of a dust storm on the daily average temperature and on the 
diurnal temperature variation, if a dust storm is in progress (TMPS 96 - TMPS123) 


If (dustM ,gt. 0.0) then 


TMPS 96 


Adjust hour of diurnal temperature maximum at 26 km 
thour 2 6 = (1. - dustA) *thour26 + dustA*16.7 
Adjust height factors for diurnal amplitudes 
Do 30 i = 1,5 

30 afact(i) = (1. - dus tA) * factor ( i ) + dustA*dfactor ( i ) 


TMPS 98 
TMPS 99 
TMPS100 
TMPS 1 0 1 
TMPS102 


Adjustments for dust storm effects on the daily average temperature for latitudes 
less than 50 degrees (absolute value) are done in TMPS 104 through TMPS110. For 
absolute latitude greater than 50 degrees, the adjustments are done in TMPS111 
through 120. 


& 


If (abslat . le . 50.)then 

dtr26 = 0.0465 + 1 . 90E-6 *abslat* *2 
T (1) = T(l) + 13 . *dustM 
T ( 2 ) = T(2) + 36 . * dustM 

T(3) = T ( 3 ) + (48. + 0 . 0024*abslat**2 ) *dustM 
T ( 4 ) = T ( 4 ) + 30 . * dustM 
T ( 5 ) = T ( 5 ) + 26 . *dustM 
Else 


alatp = 
dtr2 6 = 
4 . 394E- 
T ( 1 ) 

T { 2 ) 

T ( 3 ) 

T ( 4 ) 

T ( 5 ) 
Endif 


90. - abslat 
( <1 . 597E-6 *alatp 
3 ) *alatp 


1 . 417E-4 ) *alatp + 


T(l) + (28. - 0.3* abslat ) *dus tM 
T(2) + (11. + 0 . 015625 *alatp* *2) *dustM 
T(3) + (25. + 0 . 018125*alatp**2) * dustM 
T(4) + (50. - 0 . 4*abslat) * dustM 
T(5) + (46. - 0 . 4*abslat ) *dustM 


TMPS104 

TMPS105 

TMPS106 

TMPS107 

TMPS108 

TMPS109 

TMPS110 

TMPS111 

TMPS112 

TMPS113 

TMPS114 

TMPS115 

TMPS116 

TMPS117 

TMPS118 

TMPS119 

TMPS120 


Adjusts the diurnal amplitude from 26 to 30 km 
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ono 


dtr30 = (1. - dustA)*dtr30 + dustA*l . 072*dtr26 
Endif 


TMPS122 

TMPS123 


Steps through heights of the significant levels and evaluates the temperatures, 
T(i) including dust storm effects and diurnal variations 


Do 40 i = 1, 5 

Compute diurnal adjustment to average temperatures. 

Interpolate (extrapolate) for phase hour of maximum diurnal 
perturbation 

thour = thourO + (thour26 - thourO) *z (i) /26 . 

Make sure phase hour is in range 0 < thour <24 
If (thour .gt. 24.) thour = thour - 24. 

If (thour .It. 0.) thour = thour + 24. 

C. Longitude where maximum phase occurs 

alonmax = 15.* (thour -12.) 

C. . . Evaluate temperature including diurnal perturbation 

T(i) = T(i)*(l. + dtr3 0*afact ( i ) *cos (pil80* (dlon + alonmax))) 

40 Continue 

T(0) = Tsurf 

Re-evaluates the lapse rates (K/km) , based on the new values of temperatures at the 
significant levels 


TMPS125 

TMPS126 

TMPS127 

TMPS128 

TMPS129 

TMPS130 

TMPS131 

TMPS132 

TMPS133 

TMPS134 

TMPS135 

TMPS136 

TMPS137 

TMPS138 


Do 50 i = 1,5 

50 gam(i) = (T(i-l) - T(i))/dz(i) 


TMPS140 

TMPS141 



Subroutine : Tsur f ace 


Code : TSRF 


Description: Computes surface temperature for a given latitude and 
longitude, areocentnc longitude of the sun (Ls) and orbital radius. Also 
computes the diurnal average surface temperature, with and without dust 
storm effects. 

Called By: ATMOS 2 

Common Blocks Used (See Table 2-3 for a list of variables in coamon blocks): 

None 


Input Variables (not passed through c ommo ns ) : 


Name Type 


Description 


ALS. . 
AU. . . 
DUSTT . 
SITLA . 
SITLO . 
SUNLA . 
SUNLO . 


REAL *4 
REAL *4 
REALM 
REAL *4 
REAL *4 
REAL* 4 
REAL *4 


areocentric longitude of the sun (Ls) 

orbital radius of Mars from Sun (in astronomical units) 
relative dust storm effect (0-1) 
current latitude 
current longitude 

latitude of sub-solar point on surface 
West longitude of sub-solar point on surface 


Output Variables (not passed through commons) : 


Name Type 


Description 


TOBAR 
TAVG. 
TMAX . 
TMIN. 
TSURF 


REAL* 4 
REAL *4 
REAL *4 
REAL *4 
REAL *4 


daily average surface temperature without dust effects 
daily average surface temperature including dust effect 
daily maximum surface temperature 
daily minimum surface temperature 

local surface temperature at current position and time 


Local Variables (not passed through commons) : 


Name 

Type 

ABSLAT . . . 

REAL *4 

ABSORB . . . 

REAL *4 

ALBEDO . . . 

REAL *4 

ATRANS . . . 

REAL *4 

AVCOS . . . 

REAL *4 

A VMU .... 

REAL *4 

CAPNORTH . . 

REAL *4 

CAPSOUTH . . 

REAL *4 

CLAT .... 

REAL *4 

CLATS . . . 

REAL *4 

COSDLON . . 

REAL *4 

CRLAT . . . 

REALM 

DELTA . . . 

REAL *4 

DLO . . . . 

REAL *4 

DLON .... 

REAL *4 

DTR .... 

REALM 

F0 

REAL* 4 

PI 

REAL* 4 

PMMU .... 

REAL *4 


Description 


absolute value of latitude 

daily total solar radiation absorbed at surface 
surface reflectance for solar radiation 
factor used to compute average solar transmittance 
daily average value of the cosine of the solar zenith 
angle at the site latitude 
cosine of the solar zenith angle at solar noon 
boundary of the north polar cap 
boundary of the south polar cap 
cosine of the current latitude 
cosine of the solar latitude 

cosine of the hour-angle length of the solar day 
latitude of the boundary at which the daily total amount 
of solar radiation is 0 (total darkness for the day) 
optical depth for a vertical path through the atmosphere 
longitude difference from sun to site (also equal to the 
solar hour angle in degrees, measured from solar noon) 
hour angle length of the solar day at the site position 
n / 180 

the solar constant at current Mars orbital radius 

71 

minimum value of the cosine of the solar zenith angle 
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POLAR . 
QA. . . 


SINDLON 
SITELAT 
SLAT. . 
SLATS . 
SSALB . 
SUNLAT. 
TAMP. . 
TAUBAR. 
TLAT . . 
TLATS . 


REAL* 4 
REAL *4 


REAL*4 
REAL *4 
REAL* 4 
REAL* 4 
REAL *4 
REAL *4 
REAL*4 
REAL *4 
REAL *4 
REAL* 4 


Methodology : 


factor for temperature correction in the polar cap zone 
amplitude of the daily variation in solar heating at the 
surface (difference between average and minimum solar 
heating) 

sine of the hour-angle length of the solar day 
current latitude in radians 
sine of current latitude 
sine of the solar latitude 

single scatter albedo for atmospheric scattering 
solar latitude in radians 

amplitude of the diurnal surface temperature variation 
daily average of the solar transmittance to the surface 
tangent of the current latitude 
tangent of the solar latitude 


Computes cosine of the hour-angle length of the solar day at the site 
latitude 


cosdlon = - tlat*tlats 


TSRF 28 


Calculates crlat, the latitude at which daily insolation = 0, dlon, the length of 
the solar day, and polar, the temperature adjustment factor in the polar cap regions 
at locations where the daily insolation is zero. If | cosdlon | > 1 then either the 
insolation lasts all day (sun never sets, dlon = n radians) or there is no daylight 
(sun never rises, dlon = 0) . 


If (Abs (cosdlon) .gt. 1 . ) then 
If (sitla*sunla. It . 0 . ) Then 


TSRF 32 
TSRF 32a 


If (sunla .gt . 0 . 0)Then 
crlat = sunla - 90. 
Else 

crlat = sunla + 90 . 
Endif 


TSRF 33a 
TSRF 33b 
TSRF 33c 
TSRF 33d 
TSRF 33e 


polar = (abs (sitla-crlat) ) /25 . 19 
dlon = 0. 

Else 

dlon = pi 
Endif 

sindlon = 0. 

Else 

sindlon = sqrt(l. - cosdlon**2) 
dlon = pi/2. 

If (cosdlon. ne.0.) dlon = atan (sindlon/cosdlon) 
If (dlon .le. 0.0) dlon = pi + dlon 
Endif 


TSRF 35 
TSRF 36 
TSRF 37 
TSRF 38 
TSRF 39 
TSRF 40 
TSRF 44 
TSRF 45 
TSRF 46 
TSRF 47 
TSRF 48 
TSRF 48a 


Computes avcos, the daily average of the cosine of the solar zenith angle at the 
site latitude 


avcos = (slats* slat* dlon + clats*clat*sindlon) /pi 
If (avcos .It. 0.0) avcos = 0. 


TSRF 50 
TSRF 51 
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thS d l? St optical depth and ssalb - the single-scatter albedo to model 
values, and computes avmu, the cosine of the solar zenith angle at solar noon. 


delta =0.3 
ssalb = 0.85 
atrans = ssalb/2. 


TSRF 51b 
TSRF 51c 
TSRF 51d 


avmu = slat*slats + clat*clats 


TSRF 51f 


Computes taubar, the daily-average solar transmittance 


If (avmu . le . 0 . 0) Then 
avmu = 0.0 
taubar = 0.0 
Else 

taubar = atrans + (1. - atrans ) ♦exp ( -delta* ( 1 . /avmu) ) 
Endif 


TSRF 51h 
TSRF 51i 
TSRF 51 j 
TSRF 51k 
TSRF 511 
TSRF 51m 


Uses the Alb function to get the surface albedo 


albedo = Alb (als , sitla) 


TSRF 51o 


Calculates the daily total amount of solar radiation absorbed at the surface 


Absorb = taubar* (1. - albedo) *F0*avcos 
Finds pmmu, the cosine of the minimum of the solar zenith angle 


TSRF 53 


pmmu = slat*slats - clat*clats 
I f ( pmmu . 1 t . 0 . 0 ) pmmu = 0.0 

Evaluates Qa, the amplitude of the diurnal surface heating 

Qa = taubar* (1. - albedo) *F0* (avmu - pmmu)/2. 

Determines the limits of the northern and southern polar caps 

capnorth = 90. - polecap ( 90 . , als ) 
capsouth = -90. + polecap (-90 . , als) 


TSRF 53b 
TSRF 53c 


TSRF 53e 


TSRF 53g 
TSRF 53h 


Computes the daily average surface temperature, including the 
the polar cap region 


polar correction if in 


If (sitla. gt. capnorth .or. sitla. It. capsouth) Then 
Tavg = 140.5 + 0 . 2336 *Absorb - 8.5*polar 
Else 

Tavg = 140.5 + 0 . 8221*Absorb - 0 . 001425* (Absorb**2 ) 
Endif 


TSRF 53 j 
TSRF 53k 
TSRF 531 
TSRF 53m 
TSRF 53n 


Calculates Tamp, the amplitude of the daily surface temperature variation 
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Tamp = 0.16*Qa 


TSRF 53p 


Finds the minimum and maximum daily temperatures at site latitude 


Tmax = Tavg + Tamp 
Tmin = Tavg - Tamp 


TSRF 55 
TSRF 56 


Saves the daily average surface temperature, without dust storm effect, as TObar 


TObar = Tavg 


TSRF 59 


Compute the effects of a dust storm (if any) on Tavg, Tmax and Tmin 


If (dustT .gt. 0.0) then 
abslat = abs(sitla) 

If (abslat .le. 50.) then 

Tavg — ■ Tavg — (7. - 1 . 28E — 6*abslat**4) *dustT 
Else 

Tavg = Tavg + (7.5 - 0 . 13*abslat) *dustT 
End if 

Tamp = (0.5 - 0.25*dustT) * (Tmax - Tmin) /TObar 
Tmax = Tavg* ( 1 . + Tamp) 

Tmin = Tavg* ( 1 . - Tamp) 

Endif 


TSRF 61 
TSRF 62 
TSRF 63 
TSRF 64 
TSRF 65 
TSRF 66 
TSRF 67 
TSRF 68 
TSRF 69 
TSRF 70 
TSRF 71 


Finds the longitude difference from sun to site (equal solar hour angle, in degrees, 
measured from solar noon) , and evaluates the local surface temperature from the 
minimum and maximum daily values 


dlo = sunlo - sitlo 
if (dlo .gt. 180 . ) dlo = dlo - 360. 
if (dlo .It. -180 . ) dlo = dlo + 360. 
Tsurf = Tdiurnal (dlo, Tmin, Tmax) 


TSRF 74 
TSRF 75 
TSRF 76 
TSRF 77 
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Subroutine s Wavepert 


Code* WAVE 


Description: Zurek wave perturbations at given height, latitude and local 
time. Includes dust storm effects. Computes density perturbations from 
the original Zurek temperature wave perturbation. 

Called By: Datastep 


Common Blocks Used (See Table 2-3 for a list of variables In 


blocks) : 


WAVEDAT 


Input Variables (not passed through c ommo ns ) s 


Name Type 


CHGT . . . 
CLAT . . . 
DUSTA . . 
DUSTM . . 
1 ) 

TLOCAL . . 


REAL *4 
REAL *4 
REAL *4 
REAL* 4 

REAL *4 


Description 


current height (km) 
current latitude (degrees) 

relative magnitude of dust storm effect on phase (0-1) 
relative magnitude of dust storm effect on amplitude (0- 

local time in Mars hours (l/24th Sols) 


Output Variables (not passed through c ommo ns ) : 


Name Type 


Description 


AMP1 . . . . REAL *4 
AMP2. . . . REAL *4 
WAVE. . . . REAL * 4 


diurnal amplitude of wave (% of mean) 
semi-diurnal amplitude of wave (% of mean) 
relative density perturbation of wave (% of mean) 


Local Variables (not passed through commons): 


Name 



Type 

AMPC1 . 



REAL *4 

AMPDl . 



REAL *4 

AMPD2 . 



REAL *4 

AMPI1 . 



REAL *4 

AMPI2 . 



REAL* 4 

DLAT . . 



REAL *4 

DZ. . . 



REAL *4 

DZA . . 



REAL*4 

H . . . 



REAL *4 

11. . . 



INTEGER* 4 

12 . . . 



INTEGER* 4 

OMEGA1 . 



REAL *4 

OMTl . . 



REAL *4 

OMT2 . . 



REAL *4 

PHI . . 



REAL *4 

PH 2 . . 



REAL *4 

PHC1 . . 



REAL *4 

PHDl . . 



REAL *4 

PHD 2 . . 



REAL *4 

PHIl . . 



REAL *4 

PHI2 . . 



REAL *4 

Zl. . . 



REAL *4 

Z2 . . . 



REAL *4 


Description 


clear-atmosphere diurnal amplitude (% of mean) 

dusty-atmosphere diurnal amplitude (% of mean) 

dusty-atmosphere semi-diurnal amplitude {% of mean) 

intermediate amplitude variable used in calculations 

intermediate amplitude variable used in calculations 

relative latitude displacement from interpolation point 

intermediate value of relative height displacement 

relative height displacement from interpolation point 

height scale for vertical interpolation 

1st height index for interpolation 

2nd height index for interpolation 

diurnal frequency (15 ° / Mars hour) 

phase-adjusted diurnal frequency 

phase-adjusted semi-diurnal frequency 

diurnal phase (hours) 

semi-diurnal phase (hours) 

clear-atmosphere diurnal phase (hours) 

dusty-atmosphere diurnal phase (hours) 

dusty-atmosphere semi-diurnal phase (hours) 

intermediate phase variable used in calculations 

intermediate phase variable used in calculations 

lower height for vertical interpolation 

upper height for vertical interpolation 


Methodology : 



The Zurek wave perturbation model was adapted from Section 4.3 of ‘The Mars 
Atmosphere: Observations and Model Profiles for Mars Missions', Report No. JSC- 
24455, David E. Pitts et al., eds., and the draft report by Pitts, Tillman, Pollack 
and Zurek, 'Model Profiles of the Mars Atmosphere for the Mars Rover and Sample 
Return Mission", draft report, March 11, 1988. 

Computes diurnal frequency (15 ° / Mars hour) 


omegal = Atan ( 1 . ) / 3 . 


WAVE 13 


Evaluate il and i2, the array index values for height interpolation, the 
corresponding altitudes, zl and z2, and dza, the relative height displacement 
between zl and z2 


H = 11.26 

11 = 1 + Int ( 2 . *CHGT/H) 
If (il .It. 1) il = 1 

If (il .gt. 11) il = 11 

12 = il + 1 

zl = H* (il - 1 . ) / 2 . 
z2 = H*il/2 . 
dz = 2 . * (CHGT - zl) /H 
dza = dz 

If (dz .gt . 1 . ) dza = 1 . 


WAVE 

14 

WAVE 

15 

WAVE 

16 

WAVE 

17 

WAVE 

18 

WAVE 

19 

WAVE 

20 

WAVE 

21 

WAVE 

22 

WAVE 

23 


Uses the amplitude interpolation function (ampint) to interpolate on latitude (in 
the absolute latitude ranges 0-20°, 20-45° or > 45°) . Uses the phase interpolation 
function (phasint) to interpolate the phase. Also uses ampint and phasint to do 
height interpolation . In the interpolation section (WAVE 24 - WAVE 77), the 
amplitude variables are designated "amp' and the phase variables are designated 
«ph' , "c” denotes clear-atmosphere and "d' denotes dusty-atmosphere, 100, 120 and 

145 denote diurnal variables at latitudes 0, 20 and 45°, respectively; and 200, 220 
and 245 denote semi-diurnal variables at the same respective latitudes . 


Use the ampint function to interpolate between clear and dusty diurnal 
amplitudes, based on the value of DustM, the dust magnitude effect. Also adjust the 
semi-diurnal amplitude (applicable only to the dusty case) . 


ampl = ampint (ampcl, ampdl , DustM) /100 . 
amp 2 = Dus tM*ampd2 / 100 . 


WAVE 79 
WAVE 80 


Use the phasint function to interpolate between clear and dusty diurnal phases, 
based on the value of DustA, the magnitude of the dust effect on phase. 


phi = phasint {phcl , phd2 , DustA, omegal ) 


WAVE 81 


Compute the phase-adjusted frequencies for the diurnal and semi-diurnal components, 
and evaluate the density wave perturbation from the temperature wave model values 
[see equations (1) through (4) in the section "The Zurek Wave Perturbation Model" of 
the Release 2 Technical Report, 1993; see Appendix C] . 


omtl = omegal* (TLOCAL - phi) 
omt2 = 2 . * omegal* (TLOCAL - ph2 ) 

wave = ampl *( sin ( omtl ) - cos (omtl)) + amp2* { sin ( omt2 ) 

& cos {omt2 ) ) 


WAVE 83 
WAVE 84 
WAVE 85 
WAVE 86 
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APPENDIX B 


MARS-GRAM RELEASE #1 TECHNICAL REPORT 


This Appendix contains the technical portions of the Release #1 Report for Mars-GRAM 
(version 2.21), "The Mars Global Reference Atmospheric Model (Mars-GRAM)", Dale L. Johnson 
and Bonnie F. James (Grant Monitors), C. G. Justus and George Chimonas, October 8, 1989, prepared 
under Georgia Tech grant No. NAG8-078, for NASA Marshall Space Flight Center. Appendix 
material (giving outdated information on program input, output and running characteristics) has been 
deleted. 



ABSTRACT 


An engineering model atmosphere for Mars has been developed with many of the same 
features and capabilities of the highly successful Global Reference Atmospheric Model (GRAM) 
program for Earth's atmosphere, including mean values for density, temperature, pressure, and wind 
components, and density perturbation magnitudes and random perturbation profiles for density 
variations along specified trajectories. In the lower atmosphere of Mars (up to 75 km) the model is 
built around parameterizations of height, latitudinal, longitudinal and seasonal variations of 
temperature determined from a survey of published measurements from the Mariner and Viking 
programs. Pressure and density are inferred from the temperature by making use of the hydrostatic 
and perfect gas law relationships. For the upper atmosphere (above about 120 km), the thermospheric 
model of Stewart (1987) is used. A hydrostatic interpolation routine is used to insure a smooth 
transition from the lower portion of the model to the Stewart thermospheric model. Mars-GRAM 
includes parameterizations to simulate the effects of seasonal variation, diumal variation, dust storm 
effects, effects due to the orbital position of Mars, effects of the large seasonal variation in surface 
atmospheric pressure because of differential condensation/ sublimation of the CO 2 atmosphere in the 
polar caps, and effects of Martian atmospheric mountain wave perturbations on the magnitude of the 
expected density perturbations. The thermospheric model includes a parameterization for the effects 
of solar activity, measured by the 10.7 cm solar radio flux. Winds are computed by an aereostrophic 
(thermal wind) approximation, with the inclusion of the effects of molecular viscosity, which, because 
of the low atmospheric densities, can be very important at high altitudes. The mountain wave 
perturbation model also includes a new damping approximation due to the effects of molecular 
viscosity. 
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INTRODUCTION AND BACKGROUND 


A highly successful and well-utilized engineering model for the Earth's atmosphere, the 
Global Reference Atmospheric Model (GRAM), was developed at Georgia Tech (Justus, et al., 1975, 
1976), and has undergone several improvement cycles (Justus and Roper 1987; Justus, 1988). GRAM 
applications include orbital mechanics and lifetime studies, vehicle design and performance criteria, 
attitude control analysis problems, analysis of effects of short-term density variation from 
geomagnetic storms, and aerobraking analyses (for missions requiring return from geosynchronous 
orbit to space-station rendezvous). 

In addition to evaluating the mean density, temperature, pressure, and wind components at 
any height, latitude, longitude and monthly period, GRAM also allows for the simulation of "random 
perturbation profiles about the mean conditions. This feature permits the simulation of a large 
number of realistic density profile realizations along the same trajectory through the atmosphere, with 
realistic values of scales of variation and peak perturbation values (e.g., the random perturbation 
profiles produce values which exceed the +3 standard deviation value approximately 1 % of the time). 

With the planning activity for upcoming and proposed unmanned missions to Mars (e.g. Mars 
Observer, Mars Aeronomy Observer, Mars Rover and Sample Return), as precursors to a possible 
future manned mission, interest has developed in having a similar type of engineering-oriented 
atmospheric model as GRAM for the atmosphere of Mars. This report discusses the development of 
such a new model, the Mars Global Reference Atmospheric Model (Mars-GRAM). The Mars-GRAM 
program has been developed as a reference model atmosphere for such engineering applications as 
aerobraking or aerocapture within the atmosphere of Mars (which requires knowledge of Martian 
atmospheric density to altitudes as low as 20 km), orbiter operations at Mars (requiring atmospheric 
density information at orbital altitudes), Mars Lander entry and exit operations (requiring atmospheric 
density and wind estimates from the surface to orbital altitudes), Lander/Rover operational 
environment on the Mars surface (requiring atmospheric thermal, wind and solar radiation 
environments at the surface of Mars), and environments for design and p lannin g of balloon 
measurement systems, such as being considered on the Soviet Mars 1994 Mission (requiring 
atmospheric thermal, wind and solar radiation environments within the planetary boundary layer of 
Mars). ' ’’ 


This report outlines the observational and modeling basis behind the development of the 
Mars-GRAM program, provides documentation on the operations of the Mars-GRAM program 
(available in FORTRAN-77 on IBM-PC compatible 360k diskettes), and presents some example 
applications for the Mars-GRAM model. The report also suggests some areas for possible 
improvements for the current Mars-GRAM program, such as the incorporation of routines for the 
quantitative estimation of solar radiation at the surface of Mars, more realistic treatment of diurnal 
and latitudinal variations of density, temperature and wind at thermospheric altitudes, and the overall 
improvement of atmospheric parameterizations by further analysis of the archived data bases for the 
atmosphere of Mars. 
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THE MARS GLOBAL REFERENCE ATMOSPHERIC MODEL (MARS-GRAM) 


Mars-GRAM, is based on parameterizations to approximate, as realistically as possible, the 
tmeperature, pressure, density and winds of the Martian atmosphere, and their latitudinal, 
longitudinal, diurnal, seasonal and altitude variation, from the surface through thermospheric 
altitudes. Parameterizations are also included for the effects of global-scale dust storms on the 
variations of the thermodynamic and wind properties of the Martian atmosphere. Recently, David 
Kaplan, compiler of the definitive report "Environment of Mars, 1988" (Kaplan, 1988), has decided to 
propose Mars-GRAM as the reference model atmosphere for use by engineers on upcoming 
equipment design contracts for NASA’s Mars Rover and Sample Return Mission (Kaplan, private 
communication) . 

The near- surface air temperature on Mars is parameterized in the Mars-GRAM program by 
computing an approximate value for the geographically and seasonally-dependent daily total absorbed 
radiation flux. The daily average, maximum and minimum near-surface air temperatures are then 
calculated from a simple regression relationship assumed between the daily absorbed flux and the 
temperature parameters. The absorbed flux estimates include the variations in insulation due to the 
orbital position of Mars, the latitudinal variation of surface albedo (Pollack, et al. 198 1 ), and 
seasonally-dependent parameterizations for the polar caps (Martin and James, 1986a, 1986b; Iwasaki, 
1986; Paige and Ingersoll, 1985; Philip, 1986) and polar hood clouds (James et al., 1987). 

The resultant seasonally and latitudinally-dependent daily average, maximum and minimum 
temperatures from the Mars-GRAM model are shown in Figures 1-3. These agree well with the 
surface temperature maps of Kieffer et al. (1977), with more realistic variations in the polar regions, 
as suggested by the observations of Kieffer (1979). Seasonal variations of the daily maximum, 
minimum and average temperature at the latitudes of the Viking 1 and Viking 2 landers, as evaluated 
by Mars-GRAM, are shown in Figures 4a and 5a. These plots agree nicely with the Viking 1 and 
Viking 2 observational data reported by Ryan and Henry (1979), shown for comparison in Figures 4b 

and 5b. 


Parameterizations for the seasonal, latitudinal, diurnal and dust-storm influence on surface 
pressure were taken from data of Hess et al. (1976, 1977, and 1980), Leovy (1979), Leovy and Zurek 
(1979) Leovy (1981) and Tillman (1988). Model information on the latitudinal variations of surface 

pressure was incorporated from Haberle et al. (1982). Figure 6a presents the Mars-GRAM 

simulations for daily average surface pressure versus time at the Viking 1 and Viking 2 Lander sites. 
These compare favorably with observations of Hess et al. (1980) and Tillman (1988), as shown in 
Figure 5b (as presented by Kaplan, 1988). Parameterizations for the seasonal, latitudinal, and diurnal 
variation of pressure on the reference ellipsoid level were developed from the Viking 1 and Viking 2 
data and the model of Haberle et al. (1982). Pressures are adjusted from the values on the reference 
ellipsoid to values at the local effects of dust storms on daily mean pressure and on the amplitude ot 
the diurnal variation in pressure are evident in Figure 6b. The dust-storm effects on daily average 
pressure are assumed to be due to perturbations in the meridional circulation strength of the Martian 
Hadley cell (Haberle et al., 1982), and are therefore taken to have perturbation magnitudes which are 
antisymmetric in latitude. Thus, when a dust storm causes increased daily average surface pressure at 
northern latitudes, it is assumed to causes a comparable decrease in daily average surface pressure at 
the corresponding southern latitude. For lack of any other information, the effects of dust storms on 
the magnitude of the diurnal variation are assumed to be symmetric in latitude (i.e. equal in magnitude 
and sign at comparable northern and southern latitudes). 
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Mars-GRAM parameterizations for the geographical, seasonal and altitude dependence on 
non-dust storm, daily average temperatures above the surface come from a combination of 
observational and model values. The observational data include: Mariner 9 IRIS data from Conrath 
(1981), and Conrath data as reported by Leovy (1982) and Magalhaes (1987), Viking IRTM data 
(Martin et al„ 1982), Mariner 9 radio occultation data from Kliore et al. (1983), and Viking 1 radio 
occultation data (Fjeldbo et al., 1977), Lindal et al. (1979), and Davies (1979). Model output used 
include results from Pollack et al. (1981), Haberle et al. (1982), and Pitts et al. (1988). 

The parameterizations in Mars-GRAM for temperature variations during dust storm 
conditions were taken from observational data of Kliore et al. (1972), Jakosky and Martin (1987), 
Martin et al. (1982), Hanel et al. (1972), and Conrath (1975). Model values for dust storm effects 
were taken from Haberle et al. (1982) and Pitts et al. (1988). 

Data for the Mars-GRAM parameterizations of the height variation of the diurnal 
(longitudinal) variations of temperature about the daily mean value were taken from Viking IRTM 
observations of Martin et al. (1982), Mariner 9 IR spectroscopy data of Hanel et al. (1972), and model 
results from Zurek, as reported in Pitts et al. (1988). Diurnal temperature variations are significantly 
larger during dust storm periods than during non-dust storm conditions. 


Examples of the altitude dependence of temperature and density during non-dust storm 
conditions are provided by Figures 7 and 8, which compare Mars-GRAM simulations for the date, 
time and location of the Viking 1 Lander site with temperature observations from the Viking Lander 1 
entry profile (Seiff and Kirk, 1977) and the COSPAR Northern Hemisphere Summer Mean density 
profile (Seiff, 1982). Mars-GRAM simulations of the height and latitudinal variation of daily average 
temperatures are shown in Figure 9, for Northern Hemisphere spring equinox (areocentric longitude 
of sun, L s = 0°), in Figure 10 for Northern Hemisphere summer solstice (L s = 90°), and in Figure 1 1 
for Northern Hemisphere winter solstice (L s = 270°). Examples of Mars-GRAM simulations of dust 
storm effects on surface temperature are shown in Figures 4a and 5a. An example of a Mars-GRAM 
dust-storm height-latitude temperature cross section (at L s = 270°) is shown in Figure 12. 


In Figures 13 and 14 the model height-latitude cross sections of temperature and zonal wind 
are compared for Mars-GRAM simulation results for L s = 49° and for Mariner 9 results from L s = 43- 
54 (Leovy, 1982, from data provided by Conrath). For the same conditions as in Figure 13, Figure 15 
shows Mars-GRAM simulations of height-latitude cross sections of pressure (N/m 2 , log-base-10 scale) 
and of density (kg/m 3 , log-base- 10 scale). Figure 15a illustrates the fact that the conversion between 
height and pressure is not constant with latitude at a given seasonal time. 


The perturbations in pressure and density evident between about 30°S and 60°S in Figure 15 
are due to the influence of terrain height on surface pressure. The cross sections of Figure 15 were 
evaluated along the zero longitude meridian. Terrain height contours used in Mars-GRAM, shown in 
Figure 16, indicate that surface altitudes reach more than 3 km above the reference ellipsoid over the 
latitude range 30-60°S. The reference ellipsoid in Mars-GRAM is defined by the equatorial and polar 
radii as given by Cain et al. (1973). 


Examples of the Mars-GRAM capabilities to simulate the spatial and temporal variation of 
global dust storm effects on temperature are provided by Figures 17 and 18. These figures show both 
daily average temperatures and the range of diurnal variation in temperature during the development 
(Figure 17) and decay (Figure 18) of global dust storms. Figure 17, similar in format to Figure 5 of 
Martin et al. (1982), shows the simulated development of dust storm effects on temperature at 30 km 
altitude for a simulation of the 1977b global dust storm, which was observed by Viking IRTM. The 
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simulations in Figure 17 were based on an assumed intensity of 3.0 (the maximum possible in Mars- 
GRAM). Comparison with the data of Martin et al. indicates that better correspondence would be 
achieved by attributing an intensity of about 2.0 to the 1977b storm. Figure 18, similar in format to 
Figure 1 of Conrath (1975), is for a simulation of the decay of the 1971b global dust storm, which was 
observed by Mariner 9. The simulation in Figure 18 also assumed an intensity of 3, although the 
model would compare better with the 1971b storm observations if an intensity of about 2 were 
assumed. 

The newly-developed parameterizations for Mars-GRAM provide a simulation capability to 
altitudes reaching the base of the thermosphere. For simulations of the seasonal, geographical and 
solar-activity dependence of thermospheric conditions, the Mars-GRAM uses an adaptation of the 
thermospheric model of Stewart (Culp and Stewart, 1983, 1984; Stewart and Hanson, 1982; Stewart, 
1987; see also the revised Stewart program code given in Pitts et al., 1977). 

The Stewart model thermosphere incorporates results from a number of data and model 
sources, e.g. the Mars Reference Atmosphere (Kliore, 1982), occultation data and mass spectrometer 
data from Mariner and from Viking orbiters (Fjeldbo et al., 1966, 1970, 1977; Kliore et al., 1972; 
Stewart et al., 1972; Nier and McElroy, 1977), and data from the Viking lander atmospheric entry 
trajectories (Seiff and Kirk, 1977). It includes parameterizations to simulate the effects of solar 
activity, seasonal variation, diurnal variation magnitude, dust storm effects, and effects due to the 
orbital position of Mars. 


MARS-GRAM PROFILES OF MEAN TEMPERATURE, PRESSURE AND DENSITY 

Parameterization for the representation of surface temperature and vertical temperature 
profiles, based on summaries from observational and model data, are discussed in the previous 
section The addition of a parameterization for the variations in surface pressure, also discussed m the 
previous section, then allows for computation of vertical profiles of pressure and density by use of the 
hydrostatic and perfect gas law relations. 

Temperature profiles are built up from parameterizations of the surface temperature and 
temperature lapse rates between the significant levels set at 5, 15, 30, 50 and 75 km. Parameterization 
of the spatial and temporal variations of surface temperature were discussed in the previous section. 
From the summaries of the observational and model data, the spatial and temporal variations of the 
lapse rates between the significant levels were developed. Temperature profiles are further 
constrained to be consistent with seasonal and latitudinal variations near 25 km altitude, as presented 
in Figure 14 of Leovy (1982). 

The finally-adjusted temperature lapse rates, appropriate to the specified time and position, 
are assumed to be constant between each of the significant height levels. With the specification of 
surface pressure to serve as a boundary condition, the pressure profile is easily computed by 
integration of the hydrostatic relationship, dp/dz = -pg = -pg/RT, over each of the vertical sections of 
the constant lapse rate. After a solution for the pressure is obtained, the density is determined from 
the temperature and pressure by the perfect gas law, p = p/RT. Here R is the gas constant for a gas of 
the composition of the atmosphere of Mars and g is the altitude-dependent value of the acceleration of 
gravity on Mars (Seiff, 1982). 
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.. T k he Stewar < : model ! thermosphere starts with a specification of the temperature T F at the base of the 
ermosphere (at a height Z F ). Both T F and Z F depend on the orbital position of Mars. Between 75 km 
altitude and height Z F , the temperature profile is constructed from a combination of constant lapse rate and 
constant temperature segments as necessary to preserve the hydrostatic relationship for pressure between 
75 km and height Z F (where the Stewart model specifies the pressure to be 1.26 X 10" 4 N/m 2 ) If T P is 
greater than T(75 km), then the temperature profile between 75 km and Z F has a constant temperature [= 
T(75km>] in the lower part and a constant lapse rate in the upper rate. If T P is less than T(75 km), then the 

empera ure profile between 75 km and Z F is taken to have a constant lapse rate in the lower part and a 
constant temperature (— Tp ) in the upper part. 


TERRAIN EFFECTS ON TEMPERATURE PROFILES 


A problem in the Martian atmosphere is how to represent the effects of terrain variations on 
temperature profiles. This problem, as stated by Seiff (1982) is: 

"Near the surface, it is necessary to refer the profile to altitude above terrain as the 
terrain controls the near-surface temperatures. At higher altitudes also, it is indicated 
that radiative equilibrium is the first order control over temperature (Seiff and Kirk, 

1977) with the input at the lower boundary being radiation from the surface The 
surface temperature, in turn, is controlled by solar flux and not the surface elevation. 

Hence, we have assumed in the model that over extensive regions of depressed or 
elevated terrain ... the profile [of temperature versus height up to 100 km] . . . should 
be applied at altitudes measured above the terrain. For more local depressions or 
elevations, the model should be entered at an altitude measured relative to the mean 
surrounding surface.” 

Simply put, the Martian air on the mountain tops will not be nearly as cold compared to the air in the 
va leys as is the case on Earth. Not only is this due to the increased importance of surface radiation 
but also due to the fact that the air density is low and airflow over the mountains will not be as 
effective in moderating the effect of surface radiation on the mountaintop surface to a temperature 
near that of the approaching, high altitude, cool air flow. 

Mars-GRAM treats this problem, as does Seiff, by assuming that the surface temperature is 
unaffected by terrain variations, and that the vertical temperature profile, relative to the surface 
temperature, must be specified as a function of height above the local terrain surface. 

An example of the interrelationship between surface temperature and temperature profile is 
seen in the present Mars-GRAM results near the south pole in Figure 11. The strong latitudinal 
gradient in surface temperature near the polar cap edge is consistent with observations and other 
models. The corresponding latitudinal gradient in upper level temperatures (e g the 25 km 
temperature of Figure 14 of Leovy, on which Mars-GRAM is based) is not observed to be so large. 
Therefore there must be a strong latitudinal gradient in near-surface temperature lapse rates in this 
region near the edge of the polar cap. Details of the effects of surface radiation and terrain elevation 
on the near-surface temperature profile is an atmospheric boundary layer problem. The boundary 
layer of the Martian atmosphere has been studied only relatively little (see Sutton et al., 1978, and 
Kieffer et al., 1976). This problem of temperature effects on near surface temperatures’ and lapse rates 
is one which requires further examination in the future. 


B-7 



WIND PROFILES IN MARS-GRAM 


Wind components are evaluated in Mars-GRAM by the geostrophic (or, on Mars, 
areostrophic) approximation, obtained by a steady balance between pressure gradient and Coriolis 

forces, namely: 

u a = -(l/pf)dp/dy and v » = O/pQdP^*, (1) 

where f is the Coriolis parameter [2H sin (0), where Q. is the planetary rotation rate and 0 is the 
latitude], p is density, and the pressure gradient components are evaluated by finite differences from 
pressures evaluated at ±2.5 degrees of latitude and longitude from the position at which the wind 
components are evaluated. Since the areostrophic wind relations (1) follow the thermal wind balance, 
wind evaluated by these equations, from local horizontal pressure radients, are consistent with winds 
evaluated by vertical integration of the thermal wind equations (e.g. compare Figures 13b and 14). 

Following Fleming et al. (1988), we take the winds near the equator to be given by 

u a = -{a/[2 p ftcos(0)]}3 2 p/3y 2 , ^ 

and v a = {a/[2 p £2 cos (0)]}3 2 p/3x3y , 

where a is the radial distance from the planetary center, and the winds near the poles ( 1 0 1 > 75°) are 
computed from 

u a cos (0)/cos (75°) and v a cos (0)/cos(75°) . ( 4 ) 

Smooth transition between (1), (2) and (3) is accomplished by interpolation between 0 = 7.5° [at which 
(2) and (3) apply] and 0 = 15° [at which (1) applies]. 

At high altitudes, where the density p becomes low, the areostrophic winds become large becaus 
p is in the denominator of equations (l)-(4). To overcome this problem molecular viscosity terms are 
add ed to the momentum balance equations, of the form 

-v 3 2 u/3z 2 = -v u/L 2 and - v 3 2 v/3z 2 = -v v/L , (5) 


where v = p/p, \i is the molecular viscosity, and L is a viscous gradient scale (parameterized currently 
as being equal to the pressure scale height H). With the addition of the viscous terms in (5) to the 
balance between pressure gradient and Coriolis force terms, the viscous-modified areostrophic wind 
components become 

U V = [f 2 u a - (v/L 2 )fv a ]/[f 2 + (v/L 2 ) 2 ], ( 6a 

Vv = [f 2 v a + (v/L 2 ) f u a ] / [f 2 + (v/L 2 ) 2 ] . ( 6b 


At sufficiently high altitudes that (v/L 2 ) » f, the viscous-modified winds become 


B-8 


u v = -(L 2 /p)dp/dx and 


-(L 2 /p)a P /a y . 


Vy 


( 7 ) 


The winds evaluated by (7) do not suffer the problem of having density in the denominator, and so do 
not become too large. The winds evaluated by (7) are orthogonal to the ordinary areostrophic wind 
components, that is they are parallel to the pressure gradients instead of parallel to the isobars. This 
approximation is borne out qualitatively by Figure 2 of Bougher et al. (1988), which shows that, on a 
constant pressure surface, the winds are approximately perpendicular to the isotherms, and hence 
would be approximately parallel to the pressure gradients on a constant height surface. 


Unfortunately the parameterization of equation (6) cannot be tested thoroughly at the present, 
due to the limitations of the Stewart thermosphere model. Since Stewart’s model is intended to be 
primarily a global mean thermosphere, the horizontal gradients of pressure (as well as those of density 
and temperature) are too small to be realistic. Thus winds evaluated in the present Mars-GRAM are 
significantly smaller in magnitude than those produced by the realistic Martian thermospheric model 
of Bougher et al. (1988). 


Near the terrain surface, a simplified boundary layer model for winds is included in Mars- 
GRAM. In this model the boundary layer winds are estimated by 

u(z)= [(C D ) I/2 / k] log(z/z 0 ) u a , and (8) 

v(z)= [(C D ) I/2 / k] log(z/z D ) v a , ( 9 ) 

where the drag coefficient is estimated to be 0.0015, the von Karman constant = 0.4, and the surface 
roughness is estimated to be 0.03m. If z < 1 .6 m (the height of the Viking Lander "surface" wind 
measurements), then z is taken to be 1 .6 m in equations (8) and (9). The height at which equations (8) 
and (9) produce winds equal to u a and v a is taken to be the top of the boundary layer. 


WAVE PERTURBATION MAGNITUDES IN MARS-GRAM 

Wave-like perturbations have been observed in the Viking 1 and Viking 2 surface pressure 
data (Leovy, 1981), in the Mariner 9 IR spectroscopy data (Pirraglia and Conrath, 1974; Conrath, 

1976, 1981), and the Viking 1 and Viking 2 lander entry profiles (Seiff and Kirk, 1976, 1977). Most 
of the wave-like perturbations have been interpreted as atmospheric tides (Pirraglia and Conrath, 

1974, Conrath, 1976, Seiff, 1982), or as planetary waves (Conrath, 1981; Barnes and Hollingsworth, 
1988). The theory of thermal tides in the Martian atmosphere has been well explored (Zurek, 1976, 
1986; Zurek and Haberle, 1988). Compared to the detailed tidal analyses, relatively little attention has 
been paid to other forms of wave perturbations such as gravity waves (Gadian and Green, 1983), 
Kelvin waves and normal modes (Zurek, 1988; 1979). The very strong influence of topography on the 
wave structure in the Martian atmosphere has been documented both by observational data (Conrath, 
1976) and by theoretical analysis (Zurek, 1976). 

One of the best-known properties of gravity waves is that their amplitudes increase with 
height as the inverse square root of the mean atmospheric density. In this way, the kinetic energy of 
the wave is maintained even as the packet propagates upward into the rarer regions. The effect has 
been likened to the wave-kink that grows as it travels along a whip from the thick end near the handle 
to the thin region near the tip. The comparison is even more apt than first appears, for just as the 
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reach heights at which it becomes too large for smooth fluid behavior and it turns over into shocks 
and/or turbulence. 

One of the common sources of gravity waves is wind flow over mountains. As the 
atmospheric flow reaches the mountain, surface pressures deflect it up and over the mountain rise, and 
then restore it to level on the other side. However, the disturbance transmitted to the air flow is not 
confined to some surface region, as the pressure fields associated with the deflection have space and 
time scales that match the characteristics of gravity waves in the free atmosphere. Distorting any 
elastic medium in a way that matches its natural modes causes the excitation of waves that broadcast 
energy from the disturbance to the far reaches of the medium. So, around mountainous regions the 
atmosphere normally contains a wealth of wave activity. Mountain waves, especially the resonant 
forms known as lee waves are a familiar phenomenon of the Earth's atmosphere, and are expected to 
be of even more importance in the Martian atmosphere. 

Because erosion is much weaker in the rainless Martian environment, and because surface 
gravity is lower, surface features on Mars tend to be more severe than on Earth. Ridges and 
mountains are higher, while surface winds are comparable (and, at times, higher), so the Martian 
mountain waves provide stronger aerodynamic disturbances than are encountered on Earth. 

The waves of most interest for aerobraking or lander vehicle entry are readily modeled by the 
long- wave approximation. This component of the disturbance maps up through the atmosphere as a 
replica of the underlying terrain, and, in common with all gravity waves, the associated fractional 
density disturbances grow ever greater with altitude. Effectively the probing spacecraft will be "flying 
into the mountains", even when it is several tens of kilometers above the peaks. Limits may be 
imposed by wave saturation (Smith et al„ 1987), absorption by critical layers (Hines, 1960), or by 
dissipation of the waves by molecular viscosity (Pitteway and Hines, 1963), but over the larger terrain 
features the tops of these ghost mountains may occur as shocks or wedges of turbulence. 

In the long-wave approximation, we assume that the wave number k satisfies the inequality 


k 2 « N 2 /U 2 , (10) 

where N is the Brunt-Vaisala frequency, and U is the mean wind speed. Under this assumption the 
governing differential equation for wave perturbations is 

d 2 y/dz 2 + [N 2 /U 2 - (d 2 U/dz 2 )/U - l/4H 2 ]y(z) = 0 , (11) 

where H is the atmospheric scale height. The Fourier transform, w(z), for the vertical velocity 
perturbations is, for example, given by 

w(z)= p 0 (z)' ,/2 A(k) y(z) , (12) 

where p Q is the mean density, A(k) is presently an unknown Fourier amplitude, and y(z) is the solution 
obtained from equation (11). 

The assumption of purely tangential flow for the surface boundary condition determines the 
wave amplitude in terms of the slope dh/dx of the surface terrain and the mean wind speed U(0) that 
must be deflected over the terrain undulation. The Fourier transform of this surface condition defines 
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the component amplitudes A(k) of equation (12) through the Fourier components of the surface terrain 
shape, h(x). 

From this approach, the Fourier transform of the density variations can be obtained. Since 
there is no k dependence for the coefficients relating the Fourier transform of the density to the 
Fourier transform of the terrain surface, the Fourier inversion to obtain the vertical profile of density 
perturbations is straightforward, and yields 


Pw(*,z)/Po 


y(z)U(0)N 2 

fPoW] 

y(0)U(z)g 

lpo(z)J 


|h(x) 


(13) 


The implications of equation (13) are clear: 

(a) The horizontal structure of the density disturbance p w is an image of the terrain shape h(x). 

(b) There is the standard vertical growth factor associated with gravity waves, proportional to the 
inverse square root of the mean density. 

(c) There is a vertical wave structure y(z)/y(0), which will usually be closely sinusoidal. 

(d) There is an excitation factor [U(0)N 2 /U(z)g] relating surface conditions to the generated wave 
density field. 


These results apply equally to an isolated surface feature, such as a mountain sitting alone in 
the plains region, or an extended complex area of severe terrain. If the surface shape is known, and 
one differential equation involving the meteorological functions U(z) and N(z) is solved, the density 
perturbation, p w (x,z), at any location above the surface is obtained. Although derived here for the 
two-dimensional case, the results of equation (13) can easily be generalized to a fully three- 
dimensional terrain description. 

If non-linear processes do not limit the growth of the wave, the high kinematic molecular 
viscosity of the upper atmosphere will. Because of lower density, the molecular viscosity effect will 
be felt at lower altitudes on Mars than on Earth. The treatment of Pitteway and Hines (1963), can be 
used to model this effect. 


If a WKB approximation solution is applied for (1 1) and (13), the amplitude of the density 
perturbations due to the mountain waves can be approximated as 


|Pw( X > Z )|/p 0 


N(0)N(z) 

fPoW] 

g 

lPo(z)J 



(14) 


This form of the solution eliminates the explicit dependence on U(z) and U(0), so that all of the terms 
in (14) are readily evaluated by Mars-GRAM from the terrain height data and the thermodynamic 
parameterization models. This relation [with the terrain relief h(x) somewhat exaggerated because of 
the present coarse resolution used] is employed in Mars-GRAM to estimate the magnitude of the 
random density perturbation component below 75 km altitude. 
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In order to insure that the magnitudes estimated by (14) do not become too large [because of 
the p 0 (z) term in the denominator], a viscous damping model must also be applied here, as well as to 
the areostrophic wind components. Following the approach of Pitteway and Hines (1963), an 
approximate solution for the damping factor (appropriate in the long-wave limit), to be applied as a 
factor to equation (14), is obtained as 


expHU v N 3 H/(4 tc U 4 ) ] [1 - p 0 (z) / p o (0)]} 


( 15 ) 


where L h is the horizontal wavelength of the dominant mode of the mountain waves. L h is taken to be 
the same as the horizontal scale used in the density perturbation model (see next section). 

Relations (14) and (15) are used to estimate the density perturbation magnitudes up to a height 
of 75 km. Above the base of the thermosphere (at height Z F ), the short-term density perturbation 
magnitudes of the Stewart model are employed. Between 75 km and Z F , interpolation is used to insure 
a smooth transition between these two height regions. 

DENSITY PERTURBATION SIMULATION IN MARS-GRAM 

The same density perturbation model used in the Earth GRAM program can be utilized in 
Mars-GRAM, provided that the perturbation magnitudes, horizontal scales, and vertical scales are 
adapted to values applicable for the Martian atmosphere. The density perturbation model estimates 
relative density perturbations p'(x,y,x,t) by the relation 

p'(x+Ax,y+Ay,z+Az,t+At) = a p’(x,y,x,t) + P R(x,y,z,t) , (16) 


where a is the correlation r(Ax, Ay, Az) over the spatial separation involved between steps, given by 


r(Ax, Ay, Az) = exp{-[(Ax/Lh) 2 + (Ay/Lh) 2 + (Az/L z ) 2 ] 1/2 ) , (17) 


P is given by [1 - a 2 ] 1/2 , and R(x,y,z,t) is a normally-distributed uncorrelated random sequence, having 
the same rms magnitude as the relative density perturbations p'. 

Evaluation of the magnitudes of the density perturbations used in Mars-GRAM were 
discussed in the previous section. The vertical scale for the density perturbations, L z , was estimated 
from structure-function analysis of the Viking 1 and Viking 2 density profiles (Seiff and Kirk, 1977) 
to be 8 km, at least for altitudes up to 100 km. For the horizontal scale of the density perturbations, 
Lh, the values used in the Earth GRAM model were scaled by a factor of 1 .5, as estimated to be 
applicable for the Martian atmosphere. Namely, in Mars-GRAM we take 


L h = 30 + 0.01875 z 2 


(18) 
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with z and L h in km. An upper limit of 600 km is imposed for L* at high altitudes. 

CONCLUSIONS 

The Mars-GRAM has numerous applications as a "poor man’s global circulation model”. For 
example, the computation of all of the data necessary to describe the complete seasonal variations at 
the surface (Figures 1-3) and all altitudes (e.g., Figure 9), takes at most a few minutes on an IBM-PC 
(with 8087 co-processor; even faster on an 80286-based machine, such as an IBM-AT or compatible, 
or on an 80386-based machine). Comparable data would take many hours of computation on a 
mainframe using a 3-D global circulation model for Mars. The diurnal (longitudinal) variability 
incorporated into the Mars-GRAM program is not even available in 2-D version of a Mars global 
circulation model. 

In addition to the engineering applications envisioned for Mars-GRAM (e.g., aerocapture 
mission profile studies, Mars Rover Sample Return mission planning and design), the Mars-GRAM 
has a number of potential scientific applications. One of these is its ability to provide a realistic, 
geographically and seasonally-dependent background of temperatures and winds for studies of tides 
and the atmospheric propagation of other wave disturbances (e.g., gravity waves, mountain lee waves, 
etc.). Another application would be in providing realistic "first guess" profiles for the inversion 
processing for temperature retrievals from temperature sounders on upcoming Mars missions. 

Of course, being a parameterization model, rather than a first principles one such as a global 
circulation model, Mars-GRAM is only as good as the parameterizations built into it. Also it cannot 
test the sensitivity to variation of parameters beyond those on which the parameterizations are based 
(e.g., it cannot estimate the effects of a dust storm of twice the optical depth previously observed). 
With continued analysis of additional observational data from the Viking and Mariner programs, 
analysis of new results from global circulation models, and with new data expected to be coming in 
from the Mars Observer program, Mars-GRAM should steadily improve in its realism and reliability 
in the future. 


PLANS FOR FUTURE IMPROVEMENTS IN MARS-GRAM 

Although the goal of producing a fairly realistic model of the Martian atmosphere has been 
realized in the current version of Mars-GRAM, there is still room for substantial improvement and 
enhancement of the program. There is a considerable volume of computerized data, from the Mariner 
and Viking programs, available at government and university planetary studies laboratories around the 
country. As discussed in the previous section, since the current Mars-GRAM parameterizations were 
developed only from the limited amount of data and information in the open literature, there is room 
for considerable improvement in these parameterizations by analysis of these considerable-sized 
computer data bases. 

As discussed in the section of winds in the thermosphere, there is considerable room for 
improvement in the horizontal gradients of density, temperature and pressure in the thermospheric 
model. Simple parameterizations of the latitudinal and longitudinal (diurnal) variations of 
temperature could be developed from the realistic thermospheric simulations of the model of Bougher 
et al. (1988) (e.g., see their Figures 3 and 4). The current approach of the Stewart thermospheric 
model should then be adequate for estimating density and pressure from these improved 
thermospheric temperature estimates. 
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Since the current Mars-GRAM relies on estimates of the daily total surface solar irradiance 
for the estimation of daily maximum, average and minimum temperatures, it would be a 
straightforward addition to Mars-GRAM to add a subroutine which calculates the direct (beam) and 
diffuse (scattered) components of the time-dependent solar irradiance at the surface. Such a routine 
for estimation on the surface solar irradiance would be based on the model approach of Justus and 
Paris (1985), with parameters adapted for conditions of the Mars atmosphere and dust optical 
properties. The influence on solar irradiance due to dust storm conditions could easily be 
incorporated, since optical depth has been measured at the Viking Lander sites (Tillman et al., 1979; 
Zurek, 1982), and other dust optical properties have also been inferred (Toon et al., 1977). Addition 
of a surface solar irradiance module for Mars-GRAM would be very useful for such applications as: 
(1) analysis of thermal heating environments and solar power cell performance on surface-based 
systems, such as Mars Rover, (2) analysis of heating, life and trajectory for balloon systems, such as 
that proposed for the Soviet Mars 1994 mission, and (3) for general improvements in the capability to 
parameterize and study the Planetary Boundary Layer of Mars and its energy and momentum flux 
budgets. 

Users or potential users of Mars-GRAM who are interested in the possibility of one or more 
of these modifications being incorporated into the program are asked to contact either Dr. Justus at 
Computer Sciences Corporation or Mr. Dale Johnson at NASA Marshall Space Flight Center. 
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Figure 2 - Seasonal and latitudinal variation of daily maximum surface 
temperature, computed by the MARS-GRAM model. 
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Figure 4 - Seasonal variation of the daily maximum, mean, and 
minimum temperature at the Viking Lander 1 site (a) computed by 
MARS-GRAM, and (b) as reported by Ryan and Henry (1979) . 
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Figure 5 - Seasonal variation of the daily maximum, mean, and 
minimum temperature at the Viking Lander 2 site (a) computed 
by MARS-GRAM, and (b) as reported by Ryan and Henry (1979) . 
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Figure 6 - Seasonal variation of the daily average surface 
pressure at the Viking 1 and 2 Lander sites (a) computed 
by Mars-GRAM, and (b) as reported by Tillman (1988) . 







Figure 7 - Vertical temperature profile simulated 
by Mars-GRAM for date, time and position of Viking 
1 Lander site (solid line) and measured Viking 1 
Lander profile (dotted line). 
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Figure 8 - Vertical profile of density simulated 
by the Mars-GRAM for date, time and position of 
the Viking 1 Lander site (solid line) and the 
COSPAR model for Northern Hemisphere summer mean 
conditions (dotted line) . 
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Figure 10 - Height and latitudinal variation of daily mean 
at Northern Hemisphere summer solstice (L g = 90°) for dust- 
conditions, as computed by Mars-GRAM. 
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computed by Mars-< 








Figure 13 - Cross section of (a) temperature (K) and (b) zonal wind (m/s) 
at 9 am local time for L = 49° * 
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A3 and L = 54 (approximately early May in analogou 
n), Leovy (1982), from data provided by B. Conrath. 




Figure 15 Cross sections of (a) log-base-10 pressure (N/in ) and (b) log- 
base-10 density (kg/m ) at 9 am local time for L = 49°. 





B-36 


Figure 16 - Contours of local surface terrain height, relative to reference 
ellipsoid (km) at 10° latitude by 20° longitude resolution, 
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Figure 17 - Progression of simulated dust-storm effect on 
daily average, maximum and minimum temperature versus 
latitude and L g value (degrees) for the 1977b storm. 
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ABSTRACT 


Three major additions or modifications have been made to the Mars Global Reference 
Atmospheric Model (Mars-GRAM). (1) In addition to the interactive version, a new batch version is 
available, which uses NAMELIST input, and is completely modular, so that the main driver program 
can easily be replaced by any calling program, such as a trajectory simulation program. (2) Both the 
interactive and batch versions now have an option for treating local-scale dust storm effects, rather 
than just the global-scale dust storms in the original Mars-GRAM. (3) The Zurek wave perturbation 
model has been added, to simulate the effects of tidal perturbations, in addition to the random 
(mountain wave) perturbation model of the original Mars-GRAM. A minor modification has also 
been made which allows heights to go "below" local terrain height and return "realistic" pressure, 
density and temperature, not the surface values, as returned by the original Mars-GRAM. This feature 
will allow simulations of Mars rover paths which might go into local "valley" areas which lie below 
the average height of the present, rather coarse-resolution, terrain height data used by Mars-GRAM. 
Sample input and output of both the interactive and batch versions of Mars-GRAM are presented. 
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APPENDIX C 


MARS-GRAM RELEASE #2 TECHNICAL REPORT 


This Appendix contains the technical portions of the Release #2 Report for Mars-GRAM 
(version 3.1), "The Mars Global Reference Atmospheric Model. (Mars-GRAM) Release #2”, Bonnie 
F. James (Grant Monitor) and C. G. Justus, March 1 , 1993, prepared under Georgia Tech grant 
NAG8-877 for NASA Marshall Space Flight Center. Appendix material (giving outdated information 
on program input, output, and running characteristics) has been deleted. 



INTRODUCTION 


The Mars Global Reference Atmospheric Model (Mars-GRAM; Justus and Chimonas, 1988; 
Justus, 1990, 1991), includes simulation capabilities for mean values of density, temperature, pressure 
and wind components. Density perturbations, simulated by a mountain- wave simulation model are 
also included, for simulation of density perturbation profiles along specified trajectories through the 
atmosphere of Mars. 

Up to a height of 75 km, Mars-GRAM is based upon parameterizations of height, latitudinal, 
longitudinal and seasonal variations of temperature, determined from a survey of published 
measurements from the Mariner and Viking programs ( a complete set of references is provided by 
Justus and Chimonas, 1989). Pressure and density are inferred in Mars-GRAM from the temperature 
by making use of the hydrostatic and perfect gas law relationships. Above about 120 km, Mars- 
GRAM uses the thermospheric model of Stewart (1987). A hydrostatic interpolation routine is used 
to insure a smooth transition from the lower portion of the model to the Stewart thermosphere model. 

Mars-GRAM includes parameterizations to simulate the effects of seasonal variation, diurnal 
variation, dust storm effects, effects due to the orbital position of Mars, effects of the large seasonal 
variation in surface atmospheric pressure because of differential condensation/sublimation of the C0 2 
atmosphere in the polar caps, and effects of Martian atmospheric mountain wave perturbations on the 
magnitude of the expected density perturbations. The thermospheric model includes a 
parameterization for the effects of solar activity, measured by the 10.7 cm solar radio flux. Winds are 
computed by the thermal wind approximation, with the inclusion of the effects of molecular viscosity, 
which, because of the low atmospheric densities, can be very important at high altitudes. The 
mountain wave perturbation model also includes a damping approximation due to the effects of 
molecular viscosity. 

During much of 1990 and early 1991, the Mars Atmosphere Knowledge Requirements 
Working Group (Bourke, 1991) met by video-conference, with the objectives of 

• Understanding the sensitivity and consequences of various levels of uncertainties in the 
Martian atmosphere 

• Understanding the realistic limits on modeling the Mars atmosphere, both deterministically 
and statistically 

• Recommending a set of atmospheric information requirements to be satisfied by the robotic 
portion of the Mars Exploration Program 

Specific recommendations were made by the working group for improvements and 
additions to the Mars-GRAM program: 

(1) The addition of a capability to treat local-scale dust storms. The current Mars-GRAM treats 
all storms as growing and decaying with time, but with a size of full global dimensions. 

(2) The addition of the Zurek wave model (Pitts et al., 1988) to represent the large-scale 
temperature and density perturbations caused by atmospheric tides. The current Mars-GRAM 
perturbations are rather small-scale, gravity-wave-like variations. 
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(3) A modular version of Mars-GRAM, specifically designed for use as a subroutine in a calling 
program, such as a trajectory guidance and control analysis program. This version would 
allow easy application of user-defined perturbation models, such as sine-wave perturbations 
or hyperbolic tangent or step-like perturbations. 

The purpose of this report is to describe the Mars-GRAM modifications recently made to 
satisfy these recommendations. 


ADDITIONS TO MARS-GRAM 

The Zurek Wave Perturbation Model 

Parameters necessary to evaluate wave-structure perturbations in temperature are provided by 
Zurek’s tables on pages 1 1-12 of Pitts et al. (1988). The Zurek wave model was designed to allow 
estimation of temperature perturbations which would be produced by atmospheric tides (expected to 
be an important process in the Martian atmosphere). In order to compute wave-structure perturbations 
in density from the original temperature perturbation estimates of Zurek, a simple model is assumed, 
based on approximate hydrostatic equilibrium for these large-scale perturbations. The perfect gas 
law, p = pRT, requires that perturbations in pressure, density and temperature be related (to first 
order) by 


p’/<p> = p '/<p> + T'/<T> 


( 1 ) 


where the angle brackets denote average values. The perturbation version of the hydrostatic 
equations, 3p73z = -p'g, requires that 


3(p7<p>)3z = -(p’/<p> - p7<p>) (g/RT) = -(T7<T>) (g/RT) (2) 


If one assumes a simple cosine function for the vertical variation of T7<T>, i.e., 


T/<T> = A cos(kz) 


( 3 ) 


then equations (1) and (2) require that 

p7<p> = A [sin(kz) / (kH) - cos(kz)] 


( 4 ) 


where H is the scale height (RT/g). The temperature perturbation data provided in Zurek's table 
provide information to determine the values of k and the wave amplitude A for both clear-sky and 
dust-storm cases. If it is assumed that kH = 1, then the density perturbation may be evaluated from 
this k value, with the approximation form of equation (4), namely 



p'/<p> 


A [sin(kz) 


cos(kz)] 


(5) 


Figure 1 shows a sample evaluation of the Zurek density wave perturbation model at the 
location and time of the Viking 1 lander entry (7/20/76 12:30 GMT, at Mars latitude 22°N, longitude 
48°W). Since the random (mountain-wave) perturbations no longer constitute the only perturbations 
in Mars-GRAM, their minimum acceptable perturbation magnitudes have been decreased somewhat. 
A sample of the random density perturbations at the Viking 1 lander location and time is shown in 
Figure 2. The total density perturbations, found by adding the wave-component and random- 
component perturbations from Figures 1 and 2, are consistent with the density perturbations observed 
by the Viking 1 and 2 landers (shown as figures in the Pitts et al. 1988 report), which have peak 
values of 10-20%. 

Figure 3 shows the density wave perturbation values from the Zurek model, evaluated on a 
height-latitude cross section, through the Viking 1 longitude (48°N) at the time of the lander entry. 
Peak contour values in Figure 3 are +14% and -14% (not labeled because of the small area within 
these contours). 


Local-Scale Dust Storm Simulations 

At run time, the Mars-GRAM user selects the time of start (within seasonal bounds) for a dust 
storm (if any) to be simulated. A time profile of build-up and decay for the dust storm intensity (up to 
a selectable maximum value) is part of the program. The new additions also allow selection of a 
location (latitude and longitude) at which the dust storm is to start and a maximum radius (up to 
10,000 km) that the dust storm is allowed to grow. These parameters of dust-storm location and 
maximum radius (w) are used to compute a size factor which multiplies the intensity of the dust 
storm effects. The size factor, as a function of position and height, is given by 


size factor = 0.25 [1 + cos(90°r/r d )][l + cos(90°z/Zd)] , 


( 6 ) 


where r is the local radius from the dust storm center location, r d is the temporally varying dust-storm 
radius (up to a maximum of w), z is the local height and z d is the height of the dust storm. z d also 
grows temporally up to a maximum value of 60 km or r<j/3 (whichever is smaller). 

The dust storm radius r d and height z d are values for 1/2 the full effect. Thus the size factor 
given by equation (6) is 1 when r/r d = 0 and z/z d = 0; it is 1/2 when r/r d = 1 and z/z d = 1 ; and it is 0 
when r/r d = 2 or z/z d = 2. As illustrated by Figure 4 the function 0.5[1 + cos(x)] used in equation (6) 
is very similar to the Gaussian distribution function exp(-x 2 /Jt) frequently used as a size factor 
function is diffusion models. 

Figure 5 illustrates the effects of a local dust storm on the mean density simulated along a 
hypothetical trajectory, starting at latitude 0°N, longitude 0°W at 40 km altitude, and moving at 
constant height and latitude along longitude from 0 to 100°W, with a local dust storm centered at 
latitude 0°N, longitude 50°W, with radius = 1000 km. The values plotted in Figure 5 are the 
differences between mean density with dust storm perturbation and that with no dust storm effects. 
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Deviations in mean density of 25% or more are seen near the center position of the local storm, with 0 
deviations seen at distances from the storm center of more than twice the storm radius. 


The Mars-GRAM Batch Version 

A new batch version of Mars-GRAM has been developed, which uses NAMELIST input, and 
is completely modular, so that the main driver program can easily be replaced by any calling program, 
such as a trajectory simulation program. As with the interactive version, fixed values of trajectory 
displacements in height, latitude, longitude and time may be read in as input, or position along an 
arbitrary trajectory may be read in from an input file. In the batch version of Mars-GRAM, values of 
the trajectory displacement values can easily be changed with time by simple modifications to the 
short driver program, or from within the trajectory program which replaces the Mars-GRAM batch 
version driver program. 


ADDITIONS TO MARS-GRAM 


The wave perturbations are modeled in a new subroutine WAVEPERT and associated 
functions AMPRINT and PHASINT to interpolate wave amplitudes and phases. The Zurek wave 
model parameters are input to the program via a new BLOCK DATA routine. The local dust storm 
modifications are incorporated into the previous subroutine DUSTFACT. 

Other program modifications are that: 

(1) A minor modification has also been made which allows heights to go "below" local terrain 
height and return "realistic” pressure, density and temperature, not the surface values, as 
returned by the original Mars-GRAM. This feature will allow simulations of Mars rover 
paths which might go into local "valley" areas which lie below the average height of the 
present, rather coarse-resolution, terrain height data used by Mars-GRAM. 

(2) For simulations which are to follow the Mars-GRAM terrain heights exactly, an input height 
below -5 km will specify this option. 

(3) For 1-D plots versus either of the height variables (height above reference ellipsoid or above 
local terrain), the plotable output files have the height variable in the second (y) position 
(ordinate). This simplifies input to plot routines which do not allow run-time selection of 
which input variable is the abcissa and which is the ordinate. 


DIGITAL TERRAIN DATA 


A new set of digital terrain height data for Mars, produced by the U.S. Geological Survey, has 
been received on large format magnetic tape. The feasibility of extracting a more detailed set of 
terrain heights from this s.et, for improvement over those currently used by Mars-GRAM, has been 
examined. Some potential problems with this are that the data on DEC VAX format tapes (we have 
no easy access to a DEC VAX with 9-track tape drives), and the data on the tapes are in "pixel" form 
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for ease in producing images, not necessarily for extracting tabular data at an array (e.g., 1° X 1°) of 
fixed locations. 

This digital terrain model data base is now available on CD-ROM (Batson, 1992), and has 
been ordered in this format. However, the data were not received in time to process during this 
contract period. 
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Figure 4 - Comparison of two possible spatial distribution model functions for 
local-scale dust storms. cosf(x) - 0.5[1 + cos(x)], for -tt < x ^ with x in 
radians, and gauss (x) — exp(-x /n ) . The functions are normalized to 1 at x— 0 , 
with the area under cosf(x) from -n < x < n the same as the area under gauss (x) 


Local Dust Storm, Lon = 50 deg., Radius - 1000 
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